Revert "Store the associated R SEXP in omxMatrix instead of the localdata flag"
[openmx:openmx.git] / src / omxMatrix.c
1 /*
2  *  Copyright 2007-2013 The OpenMx Project
3  *
4  *  Licensed under the Apache License, Version 2.0 (the "License");
5  *  you may not use this file except in compliance with the License.
6  *  You may obtain a copy of the License at
7  *
8  *       http://www.apache.org/licenses/LICENSE-2.0
9  *
10  *  Unless required by applicable law or agreed to in writing, software
11  *  distributed under the License is distributed on an "AS IS" BASIS,
12  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  *  See the License for the specific language governing permissions and
14  *  limitations under the License.
15  */
16
17 /***********************************************************
18 *
19 *  omxMatrix.cc
20 *
21 *  Created: Timothy R. Brick    Date: 2008-11-13 12:33:06
22 *
23 *       Contains code for the omxMatrix class
24 *   omxDataMatrices hold necessary information to simplify
25 *       dealings between the OpenMX back end and BLAS.
26 *
27 **********************************************************/
28 #include "omxMatrix.h"
29
30 // forward declarations
31 omxMatrix* fillMatrixHelperFunction(omxMatrix* om, SEXP matrix, omxState* state,
32         unsigned short hasMatrixNumber, int matrixNumber);
33
34 const char omxMatrixMajorityList[3] = "Tn";             // BLAS Column Majority.
35
36 void omxPrintMatrix(omxMatrix *source, char* header) {
37         int j, k;
38
39         Rprintf("%s: (%d x %d) [%s-major]\n", header, source->rows, source->cols, (source->colMajor?"col":"row"));
40         if(OMX_DEBUG_MATRIX) {Rprintf("Matrix Printing is at %0x\n", source);}
41
42         if(source->colMajor) {
43                 for(j = 0; j < source->rows; j++) {
44                         for(k = 0; k < source->cols; k++) {
45                                 Rprintf("\t%3.6f", source->data[k*source->rows+j]);
46                         }
47                         Rprintf("\n");
48                 }
49         } else {
50                 for(j = 0; j < source->cols; j++) {
51                         for(k = 0; k < source->rows; k++) {
52                                 Rprintf("\t%3.6f", source->data[k*source->cols+j]);
53                         }
54                         Rprintf("\n");
55                 }
56         }
57 }
58
59 omxMatrix* omxInitMatrix(omxMatrix* om, int nrows, int ncols, unsigned short isColMajor, omxState* os) {
60
61         if(om == NULL) om = (omxMatrix*) R_alloc(1, sizeof(omxMatrix));
62         if(OMX_DEBUG_MATRIX) { Rprintf("Initializing matrix 0x%0x to (%d, %d) with state at 0x%x.\n", om, nrows, ncols, os); }
63
64         om->hasMatrixNumber = 0;
65         om->rows = nrows;
66         om->cols = ncols;
67         om->colMajor = (isColMajor ? 1 : 0);
68
69         om->originalRows = om->rows;
70         om->originalCols = om->cols;
71         om->originalColMajor=om->colMajor;
72
73         if(om->rows == 0 || om->cols == 0) {
74                 om->data = NULL;
75                 om->localData = FALSE;
76         } else {
77                 om->data = (double*) Calloc(nrows * ncols, double);
78                 om->localData = TRUE;
79         }
80
81         om->populateFrom = NULL;
82         om->populateFromCol = NULL;
83         om->populateFromRow = NULL;
84         om->populateToCol = NULL;
85         om->populateToRow = NULL;
86
87         om->numPopulateLocations = 0;
88
89         om->aliasedPtr = NULL;
90         om->algebra = NULL;
91         om->fitFunction = NULL;
92
93         om->currentState = os;
94         om->isTemporary = FALSE;
95         om->name = NULL;
96         om->isDirty = FALSE;
97
98         omxMatrixLeadingLagging(om);
99
100         return om;
101
102 }
103
104 omxMatrix* omxInitTemporaryMatrix(omxMatrix* om, int nrows, int ncols, unsigned short isColMajor, omxState* os) {
105
106         if(om == NULL) {
107                 om = (omxMatrix*) Calloc(1, omxMatrix);
108         }
109
110         om = omxInitMatrix(om, nrows, ncols, isColMajor, os);
111         om->isTemporary = TRUE;
112         
113         return(om);
114
115 }
116
117 void omxCopyMatrix(omxMatrix *dest, omxMatrix *orig) {
118         /* Copy a matrix.  NOTE: Matrix maintains its algebra bindings. */
119
120         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("omxCopyMatrix"); }
121
122         int regenerateMemory = TRUE;
123         int numPopLocs = orig->numPopulateLocations;
124
125         if(dest->localData && (dest->originalRows == orig->rows && dest->originalCols == orig->cols)) {
126                 regenerateMemory = FALSE;                               // If it's local data and the right size, we can keep memory.
127         }
128
129         dest->rows = orig->rows;
130         dest->cols = orig->cols;
131         dest->colMajor = orig->colMajor;
132         dest->originalRows = dest->rows;
133         dest->originalCols = dest->cols;
134         dest->originalColMajor = dest->colMajor;
135         dest->isDirty = orig->isDirty;
136
137         dest->numPopulateLocations = numPopLocs;
138         if (numPopLocs > 0) {
139                 dest->populateFrom = (int*)R_alloc(numPopLocs, sizeof(int));
140                 dest->populateFromRow = (int*)R_alloc(numPopLocs, sizeof(int));
141                 dest->populateFromCol = (int*)R_alloc(numPopLocs, sizeof(int));
142                 dest->populateToRow = (int*)R_alloc(numPopLocs, sizeof(int));
143                 dest->populateToCol = (int*)R_alloc(numPopLocs, sizeof(int));
144                 
145                 memcpy(dest->populateFrom, orig->populateFrom, numPopLocs * sizeof(int));
146                 memcpy(dest->populateFromRow, orig->populateFromRow, numPopLocs * sizeof(int));
147                 memcpy(dest->populateFromCol, orig->populateFromCol, numPopLocs * sizeof(int));
148                 memcpy(dest->populateToRow, orig->populateToRow, numPopLocs * sizeof(int));
149                 memcpy(dest->populateToCol, orig->populateToCol, numPopLocs * sizeof(int));
150         }
151
152         if(dest->rows == 0 || dest->cols == 0) {
153                 omxFreeMatrixData(dest);
154                 dest->data = NULL;
155                 dest->localData=FALSE;
156         } else {
157                 if(regenerateMemory) {
158                         omxFreeMatrixData(dest);                                                                                        // Free and regenerate memory
159                         dest->data = (double*) Calloc(dest->rows * dest->cols, double);
160                 }
161                 memcpy(dest->data, orig->data, dest->rows * dest->cols * sizeof(double));
162                 dest->localData = TRUE;
163         }
164
165         dest->aliasedPtr = NULL;
166
167         omxMatrixLeadingLagging(dest);
168 }
169
170 void omxAliasMatrix(omxMatrix *dest, omxMatrix *src) {
171         omxCopyMatrix(dest, src);
172         dest->aliasedPtr = src;                                 // Alias now follows back matrix precisely.
173 }
174
175 void omxFreeMatrixData(omxMatrix * om) {
176
177         if(om->localData && om->data != NULL) {
178                 if(OMX_DEBUG_MATRIX) { Rprintf("Freeing matrix at 0x%0x. Localdata = %d.\n", om->data, om->localData); }
179                 Free(om->data);
180                 om->data = NULL;
181                 om->localData = FALSE;
182         }
183
184 }
185
186 void omxFreeAllMatrixData(omxMatrix *om) {
187     
188     if(om == NULL) return;
189
190         if(OMX_DEBUG) { 
191             Rprintf("Freeing matrix at 0x%0x with data = 0x%x, algebra 0x%x, and fit function 0x%x.\n", 
192                 om, om->data, om->algebra); 
193         }
194
195         if(om->localData && om->data != NULL) {
196                 Free(om->data);
197                 om->data = NULL;
198                 om->localData = FALSE;
199         }
200
201         if(om->algebra != NULL) {
202                 omxFreeAlgebraArgs(om->algebra);
203                 om->algebra = NULL;
204         }
205
206         if(om->fitFunction != NULL) {
207                 omxFreeFitFunctionArgs(om->fitFunction);
208                 om->fitFunction = NULL;
209         }
210         
211         if(om->isTemporary) {
212                 Free(om);
213                 om = NULL;
214         }
215
216 }
217
218 /**
219  * Copies an omxMatrix to a new R matrix object
220  *
221  * \param om the omxMatrix to copy
222  * \return a PROTECT'd SEXP for the R matrix object
223  */
224 SEXP omxExportMatrix(omxMatrix *om) {
225         SEXP nextMat;
226         PROTECT(nextMat = allocMatrix(REALSXP, om->rows, om->cols));
227         for(int row = 0; row < om->rows; row++) {
228                 for(int col = 0; col < om->cols; col++) {
229                         REAL(nextMat)[col * om->rows + row] =
230                                 omxMatrixElement(om, row, col);
231                 }
232         }
233         return nextMat;
234 }
235
236 void omxZeroByZeroMatrix(omxMatrix *om) {
237         if (om->rows > 0 || om->cols > 0) {
238                 omxResizeMatrix(om, 0, 0, FALSE);
239         }
240 }
241
242 omxMatrix* omxNewIdentityMatrix(int nrows, omxState* state) {
243         omxMatrix* newMat = NULL;
244         int l,k;
245
246         newMat = omxInitMatrix(newMat, nrows, nrows, FALSE, state);
247         for(k =0; k < newMat->rows; k++) {
248                 for(l = 0; l < newMat->cols; l++) {
249                         if(l == k) {
250                                 omxSetMatrixElement(newMat, k, l, 1);
251                         } else {
252                                 omxSetMatrixElement(newMat, k, l, 0);
253                         }
254                 }
255         }
256         return newMat;
257 }
258
259 omxMatrix* omxDuplicateMatrix(omxMatrix* src, omxState* newState) {
260         omxMatrix* newMat;
261     
262         if(src == NULL) return NULL;
263         newMat = omxInitMatrix(NULL, src->rows, src->cols, FALSE, newState);
264         omxCopyMatrix(newMat, src);
265         newMat->hasMatrixNumber = src->hasMatrixNumber;
266         newMat->matrixNumber    = src->matrixNumber;
267     
268     return newMat;    
269 }
270
271 void omxResizeMatrix(omxMatrix *om, int nrows, int ncols, unsigned short keepMemory) {
272         // Always Recompute() before you Resize().
273         if(OMX_DEBUG_MATRIX) { 
274                 Rprintf("Resizing matrix from (%d, %d) to (%d, %d) (keepMemory: %d)", 
275                         om->rows, om->cols, 
276                         nrows, ncols, keepMemory);
277         }
278         if((keepMemory == FALSE) && (om->rows != nrows || om->cols != ncols)) {
279                 if(OMX_DEBUG_MATRIX) { Rprintf(" and regenerating memory to do it"); }
280                 omxFreeMatrixData(om);
281                 om->data = (double*) Calloc(nrows * ncols, double);
282                 om->localData = TRUE;
283         } else if(om->originalRows * om->originalCols < nrows * ncols) {
284                 warning("Upsizing an existing matrix may cause undefined behavior.\n"); // TODO: Define this behavior?
285         }
286
287         if(OMX_DEBUG_MATRIX) { Rprintf(".\n"); }
288         om->rows = nrows;
289         om->cols = ncols;
290         if(keepMemory == FALSE) {
291                 om->originalRows = om->rows;
292                 om->originalCols = om->cols;
293         }
294
295         omxMatrixLeadingLagging(om);
296 }
297
298 void omxResetAliasedMatrix(omxMatrix *om) {
299         om->rows = om->originalRows;
300         om->cols = om->originalCols;
301         if(om->aliasedPtr != NULL) {
302                 memcpy(om->data, om->aliasedPtr->data, om->rows*om->cols*sizeof(double));
303                 om->colMajor = om->aliasedPtr->colMajor;
304         }
305         omxMatrixLeadingLagging(om);
306 }
307
308 double* omxLocationOfMatrixElement(omxMatrix *om, int row, int col) {
309         int index = 0;
310         if(om->colMajor) {
311                 index = col * om->rows + row;
312         } else {
313                 index = row * om->cols + col;
314         }
315         return om->data + index;
316 }
317
318 void vectorElementError(int index, int numrow, int numcol) {
319         char *errstr = calloc(250, sizeof(char));
320         if ((numrow > 1) && (numcol > 1)) {
321                 sprintf(errstr, "Requested improper index (%d) from a malformed vector of dimensions (%d, %d).", 
322                         index, numrow, numcol);
323         } else {
324                 int length = (numrow > 1) ? numrow : numcol;
325                 sprintf(errstr, "Requested improper index (%d) from vector of length (%d).", 
326                         index, length);
327         }
328         error(errstr);
329         free(errstr);  // TODO not reached
330 }
331
332 void setMatrixError(omxMatrix *om, int row, int col, int numrow, int numcol) {
333         char *errstr = calloc(250, sizeof(char));
334         static const char *matrixString = "matrix";
335         static const char *algebraString = "algebra";
336         static const char *fitString = "fit function";
337         const char *typeString;
338         if (om->algebra != NULL) {
339                 typeString = algebraString;
340         } else if (om->fitFunction != NULL) {
341                 typeString = fitString;
342         } else {
343                 typeString = matrixString;
344         }
345         if (om->name == NULL) {
346                 sprintf(errstr, "Attempted to set row and column (%d, %d) in %s with dimensions %d x %d.", 
347                         row, col, typeString, numrow, numcol);
348         } else {
349                 sprintf(errstr, "Attempted to set row and column (%d, %d) in %s \"%s\" with dimensions %d x %d.", 
350                         row, col, typeString, om->name, numrow, numcol);
351         }
352         error(errstr);
353         free(errstr);  // TODO not reached
354 }
355
356 void matrixElementError(int row, int col, int numrow, int numcol) {
357         char *errstr = calloc(250, sizeof(char));
358         sprintf(errstr, "Requested improper value (%d, %d) from (%d, %d) matrix.",
359                 row, col, numrow, numcol);
360         error(errstr);
361         free(errstr);  // TODO not reached
362 }
363
364 void setVectorError(int index, int numrow, int numcol) {
365         char *errstr = calloc(250, sizeof(char));
366         if ((numrow > 1) && (numcol > 1)) {
367                 sprintf(errstr, "Attempting to set improper index (%d) from a malformed vector of dimensions (%d, %d).", 
368                         index, numrow, numcol);
369         } else {
370                 int length = (numrow > 1) ? numrow : numcol;
371                 sprintf(errstr, "Setting improper index (%d) from vector of length %d.", 
372                         index, length);
373         }
374         error(errstr);
375         free(errstr);  // TODO not reached
376 }
377
378 double omxAliasedMatrixElement(omxMatrix *om, int row, int col) {
379         int index = 0;
380         if(row >= om->originalRows || col >= om->originalCols) {
381                 char *errstr = calloc(250, sizeof(char));
382                 sprintf(errstr, "Requested improper value (%d, %d) from (%d, %d) matrix.", 
383                         row + 1, col + 1, om->originalRows, om->originalCols);
384                 error(errstr);
385                 free(errstr);  // TODO not reached
386         return (NA_REAL);
387         }
388         if(om->colMajor) {
389                 index = col * om->originalRows + row;
390         } else {
391                 index = row * om->originalCols + col;
392         }
393         return om->data[index];
394 }
395
396 void omxMarkDirty(omxMatrix *om) { om->isDirty = TRUE; }
397 void omxMarkClean(omxMatrix *om) { om->isDirty = FALSE; }
398
399 unsigned short omxMatrixNeedsUpdate(omxMatrix *om) {
400         for(int i = 0; i < om->numPopulateLocations; i++) {
401                 int index = om->populateFrom[i];
402                 omxMatrix* sourceMatrix;
403                 if (index < 0) {
404                         sourceMatrix = om->currentState->matrixList[~index];
405                 } else {
406                         sourceMatrix = om->currentState->algebraList[index];
407                 }
408                 if(omxNeedsUpdate(sourceMatrix)) return TRUE;   // Make sure it's up to date
409         }
410     return FALSE;
411 };
412
413 omxMatrix* omxNewMatrixFromRPrimitive(SEXP rObject, omxState* state, 
414         unsigned short hasMatrixNumber, int matrixNumber) {
415 /* Creates and populates an omxMatrix with details from an R matrix object. */
416         omxMatrix *om = NULL;
417         om = omxInitMatrix(NULL, 0, 0, FALSE, state);
418         return omxFillMatrixFromRPrimitive(om, rObject, state, hasMatrixNumber, matrixNumber);
419 }
420
421 omxMatrix* omxFillMatrixFromRPrimitive(omxMatrix* om, SEXP rObject, omxState* state,
422         unsigned short hasMatrixNumber, int matrixNumber) {
423 /* Populates the fields of a omxMatrix with details from an R object. */
424         if(!isMatrix(rObject) && !isVector(rObject)) { // Sanity Check
425                 error("Recieved unknown matrix type in omxFillMatrixFromRPrimitive.");
426         }
427         return(fillMatrixHelperFunction(om, rObject, state, hasMatrixNumber, matrixNumber));
428 }
429
430
431
432 omxMatrix* fillMatrixHelperFunction(omxMatrix* om, SEXP matrix, omxState* state,
433         unsigned short hasMatrixNumber, int matrixNumber) {
434
435         SEXP matrixDims;
436         int* dimList;
437
438         if(OMX_DEBUG) { Rprintf("Filling omxMatrix from R matrix.\n"); }
439
440         if(om == NULL) {
441                 om = omxInitMatrix(NULL, 0, 0, FALSE, state);
442         }
443
444         om->data = REAL(AS_NUMERIC(matrix));    // TODO: Class-check first?
445
446         if(isMatrix(matrix)) {
447                 PROTECT(matrixDims = getAttrib(matrix, R_DimSymbol));
448                 dimList = INTEGER(matrixDims);
449                 om->rows = dimList[0];
450                 om->cols = dimList[1];
451                 UNPROTECT(1);   // MatrixDims
452         } else if (isVector(matrix)) {          // If it's a vector, assume it's a row vector. BLAS doesn't care.
453                 if(OMX_DEBUG) { Rprintf("Vector discovered.  Assuming rowity.\n"); }
454                 om->rows = 1;
455                 om->cols = length(matrix);
456         }
457         if(OMX_DEBUG) { Rprintf("Matrix connected to (%d, %d) matrix or MxMatrix.\n", om->rows, om->cols); }
458
459         om->localData = FALSE;
460         om->colMajor = TRUE;
461         om->originalRows = om->rows;
462         om->originalCols = om->cols;
463         om->originalColMajor = TRUE;
464         om->aliasedPtr = NULL;
465         om->algebra = NULL;
466         om->fitFunction = NULL;
467         om->currentState = state;
468         om->hasMatrixNumber = hasMatrixNumber;
469         om->matrixNumber = matrixNumber;
470         om->isDirty = FALSE;
471
472         if(OMX_DEBUG) { Rprintf("Pre-compute call.\n");}
473         omxMatrixLeadingLagging(om);
474         if(OMX_DEBUG) { Rprintf("Post-compute call.\n");}
475
476         if(OMX_DEBUG) {
477                 omxPrintMatrix(om, "Finished importing matrix");
478         }
479
480         return om;
481 }
482
483 void omxProcessMatrixPopulationList(omxMatrix* matrix, SEXP matStruct) {
484
485         if(OMX_DEBUG) { Rprintf("Processing Population List: %d elements.\n", length(matStruct) - 1); }
486         SEXP subList;
487
488         if(length(matStruct) > 1) {
489                 int numPopLocs = length(matStruct) - 1;
490                 matrix->numPopulateLocations = numPopLocs;
491                 matrix->populateFrom = (int*)R_alloc(numPopLocs, sizeof(int));
492                 matrix->populateFromRow = (int*)R_alloc(numPopLocs, sizeof(int));
493                 matrix->populateFromCol = (int*)R_alloc(numPopLocs, sizeof(int));
494                 matrix->populateToRow = (int*)R_alloc(numPopLocs, sizeof(int));
495                 matrix->populateToCol = (int*)R_alloc(numPopLocs, sizeof(int));
496         }
497
498         for(int i = 0; i < length(matStruct)-1; i++) {
499                 PROTECT(subList = AS_INTEGER(VECTOR_ELT(matStruct, i+1)));
500
501                 int* locations = INTEGER(subList);
502                 if(OMX_DEBUG) { Rprintf("."); } //:::
503                 matrix->populateFrom[i] = locations[0];
504                 matrix->populateFromRow[i] = locations[1];
505                 matrix->populateFromCol[i] = locations[2];
506                 matrix->populateToRow[i] = locations[3];
507                 matrix->populateToCol[i] = locations[4];
508
509                 UNPROTECT(1); // subList
510         }
511 }
512
513 void omxToggleRowColumnMajor(omxMatrix *mat) {
514
515         int i, j;
516         int nrows = mat->rows;
517         int ncols = mat->cols;
518         
519         double *newdata = (double*) Calloc(nrows * ncols, double);
520         double *olddata = mat->data;
521
522         if (mat->colMajor) {
523                 for(i = 0; i < ncols; i++) {
524                         for(j = 0; j < nrows; j++) {
525                                 newdata[i + ncols * j] = olddata[i * nrows + j];
526                         }
527                 }
528         } else {
529                 for(i = 0; i < nrows; i++) {
530                         for(j = 0; j < ncols; j++) {
531                                 newdata[i + nrows * j] = olddata[i * ncols + j];
532                         }
533                 }
534         }
535
536         if (mat->localData) {
537                 Free(mat->data);
538         }
539
540         mat->localData = TRUE;
541         mat->data = newdata;
542         mat->colMajor = !mat->colMajor;
543 }
544
545 void omxTransposeMatrix(omxMatrix *mat) {
546         mat->colMajor = !mat->colMajor;
547         
548         if(mat->rows != mat->cols){
549         int mid = mat->rows;
550         mat->rows = mat->cols;
551         mat->cols = mid;
552         }
553         
554         omxMatrixLeadingLagging(mat);
555 }
556
557 void omxRemoveElements(omxMatrix *om, int numRemoved, int removed[]) {
558
559         if(numRemoved < 1) { return; }
560
561         int oldElements;
562
563         if (om->rows > 1) {
564                 if(om->aliasedPtr == NULL) {
565                         if(om->originalRows == 0) {
566                                 om->originalRows = om->rows;
567                         }
568                         oldElements = om->originalRows;
569                 } else {
570                         oldElements = om->aliasedPtr->rows;
571                 }
572                 om->rows = oldElements - numRemoved;
573         } else {
574                 if(om->aliasedPtr == NULL) {
575                         if(om->originalCols == 0) {
576                                 om->originalCols = om->cols;
577                         }
578                         oldElements = om->originalCols;
579                 } else {
580                         oldElements = om->aliasedPtr->cols;
581                 }
582                 om->cols = oldElements - numRemoved;
583         }
584
585         int nextElement = 0;
586
587         for(int j = 0; j < oldElements; j++) {
588                 if(!removed[j]) {
589                         if(om->aliasedPtr == NULL) {
590                                 omxUnsafeSetVectorElement(om, nextElement, omxUnsafeVectorElement(om, j));
591                         } else {
592                                 omxUnsafeSetVectorElement(om, nextElement, omxUnsafeVectorElement(om->aliasedPtr, j));
593                         }
594                         nextElement++;
595                 }
596         }
597
598         omxMatrixLeadingLagging(om);
599 }
600
601 void omxRemoveRowsAndColumns(omxMatrix *om, int numRowsRemoved, int numColsRemoved, int rowsRemoved[], int colsRemoved[])
602 {
603     // TODO: Create short-circuit form of omxRemoveRowsAndCols to remove just rows or just columns.
604 //      if(OMX_DEBUG_MATRIX) { Rprintf("Removing %d rows and %d columns from 0x%0x.\n", numRowsRemoved, numColsRemoved, om);}
605
606         if(numRowsRemoved < 1 && numColsRemoved < 1) { return; }
607
608         int oldRows, oldCols;
609
610         if(om->aliasedPtr == NULL) {
611                 if(om->originalRows == 0 || om->originalCols == 0) {
612                         om->originalRows = om->rows;
613                         om->originalCols = om->cols;
614                 }
615                 oldRows = om->originalRows;
616                 oldCols = om->originalCols;
617         } else {
618                 oldRows = om->aliasedPtr->rows;
619                 oldCols = om->aliasedPtr->cols;
620         }
621
622         int nextCol = 0;
623         int nextRow = 0;
624
625         if(om->rows > om->originalRows || om->cols > om->originalCols) {        // sanity check.
626                 error("Aliased Matrix is too small for alias.");
627         }
628
629         om->rows = oldRows - numRowsRemoved;
630         om->cols = oldCols - numColsRemoved;
631
632         // Note:  This really aught to be done using a matrix multiply.  Why isn't it?
633         for(int j = 0; j < oldCols; j++) {
634                 if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Handling column %d/%d...", j, oldCols);}
635                 if(colsRemoved[j]) {
636                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Removed.\n");}
637                         continue;
638                 } else {
639                         nextRow = 0;
640                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Rows (max %d): ", oldRows); }
641                         for(int k = 0; k < oldRows; k++) {
642                                 if(rowsRemoved[k]) {
643                                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("%d removed....", k);}
644                                         continue;
645                                 } else {
646                                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("%d kept....", k);}
647                                         if(om->aliasedPtr == NULL) {
648                                                 if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Self-aliased matrix access.\n");}
649                                                 omxSetMatrixElement(om, nextRow, nextCol, omxAliasedMatrixElement(om, k, j));
650                                         } else {
651                                                 if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Matrix 0x%x re-aliasing to 0x%x.\n", om, om->aliasedPtr);}
652                                                 omxSetMatrixElement(om, nextRow, nextCol, omxMatrixElement(om->aliasedPtr, k,  j));
653                                         }
654                                         nextRow++;
655                                 }
656                         }
657                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("\n");}
658                         nextCol++;
659                 }
660         }
661
662         omxMatrixLeadingLagging(om);
663 }
664
665 /* Function wrappers that switch based on inclusion of algebras */
666 void omxPrint(omxMatrix *source, char* d) {                                     // Pretty-print a (small) matrix
667     if(source == NULL) Rprintf("%s is NULL.\n", d);
668         else if(source->algebra != NULL) omxAlgebraPrint(source->algebra, d);
669         else if(source->fitFunction != NULL) omxFitFunctionPrint(source->fitFunction, d);
670         else omxPrintMatrix(source, d);
671 }
672
673 void omxPopulateSubstitutions(omxMatrix *om) {
674         for(int i = 0; i < om->numPopulateLocations; i++) {
675                 int index = om->populateFrom[i];
676                 omxMatrix* sourceMatrix;
677                 if (index < 0) {
678                         sourceMatrix = om->currentState->matrixList[~index];
679                 } else {
680                         sourceMatrix = om->currentState->algebraList[index];
681                 }
682                 if (sourceMatrix != NULL) {
683                         omxRecompute(sourceMatrix);                             // Make sure it's up to date
684                         double value = omxMatrixElement(sourceMatrix, om->populateFromRow[i], om->populateFromCol[i]);
685                         omxSetMatrixElement(om, om->populateToRow[i], om->populateToCol[i], value);
686                 }
687         }
688 }
689
690 void omxMatrixLeadingLagging(omxMatrix *om) {
691         om->majority = &(omxMatrixMajorityList[(om->colMajor?1:0)]);
692         om->minority = &(omxMatrixMajorityList[(om->colMajor?0:1)]);
693         om->leading = (om->colMajor?om->rows:om->cols);
694         om->lagging = (om->colMajor?om->cols:om->rows);
695 }
696
697 unsigned short omxNeedsUpdate(omxMatrix *matrix) {
698         if (OMX_DEBUG_ALGEBRA) 
699                 Rprintf("Matrix 0x%0x needs update: ", matrix);
700         if (matrix->hasMatrixNumber && !matrix->isDirty) {
701                 if (OMX_DEBUG_ALGEBRA) Rprintf("No\n");
702                 return(FALSE);
703         } else {
704                 if (OMX_DEBUG_ALGEBRA) Rprintf("Yes\n");
705                 return(TRUE);
706         }
707 }
708
709 void omxRecompute(omxMatrix *matrix) {
710         if(matrix->numPopulateLocations > 0) omxPopulateSubstitutions(matrix);
711         else if(!omxNeedsUpdate(matrix)) /* do nothing */;
712         else if(matrix->algebra != NULL) omxAlgebraRecompute(matrix->algebra);
713         else if(matrix->fitFunction != NULL) {
714                 omxFitFunctionCompute(matrix->fitFunction, FF_COMPUTE_FIT, NULL);
715         }
716 }
717
718 void omxCompute(omxMatrix *matrix) {
719         if(matrix->numPopulateLocations > 0) omxPopulateSubstitutions(matrix);
720         else if (matrix->algebra != NULL) omxAlgebraCompute(matrix->algebra);
721         else if(matrix->fitFunction != NULL) {
722                 omxFitFunctionCompute(matrix->fitFunction, FF_COMPUTE_FIT, NULL);
723         }
724 }