Committing fixes to multicore execution. Added argument 'copyState' to omxCopyMatrix...
[openmx:openmx.git] / src / omxMatrix.c
1 /*
2  *  Copyright 2007-2009 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->rows = nrows;
65         om->cols = ncols;
66         om->colMajor = (isColMajor?1:0);
67
68         om->originalRows = om->rows;
69         om->originalCols = om->cols;
70         om->originalColMajor=om->colMajor;
71
72         if(om->rows == 0 || om->cols == 0) {
73                 om->data = NULL;
74                 om->localData = FALSE;
75         } else {
76                 om->data = (double*) Calloc(nrows * ncols, double);
77                 om->localData = TRUE;
78         }
79
80         om->populateFrom = NULL;
81         om->populateFromCol = NULL;
82         om->populateFromRow = NULL;
83         om->populateToCol = NULL;
84         om->populateToRow = NULL;
85
86         om->numPopulateLocations = 0;
87
88         om->aliasedPtr = NULL;
89         om->algebra = NULL;
90         om->objective = NULL;
91
92         om->currentState = os;
93         om->lastCompute = -2;
94         om->lastRow = -2;
95         om->isTemporary = FALSE;
96
97         omxMatrixCompute(om);
98
99         return om;
100
101 }
102
103 omxMatrix* omxInitTemporaryMatrix(omxMatrix* om, int nrows, int ncols, unsigned short isColMajor, omxState* os) {
104
105         if(om == NULL) {
106                 om = (omxMatrix*) Calloc(1, omxMatrix);
107         }
108
109         om = omxInitMatrix(om, nrows, ncols, isColMajor, os);
110         om->isTemporary = TRUE;
111         
112         return(om);
113
114 }
115
116 void omxCopyMatrix(omxMatrix *dest, omxMatrix *orig, int copyState) {
117         /* Duplicate a matrix.  NOTE: Matrix maintains its algebra bindings. */
118
119         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("omxCopyMatrix"); }
120
121         int regenerateMemory = TRUE;
122
123         if(dest->localData && (dest->originalRows == orig->rows && dest->originalCols == orig->cols)) {
124                 regenerateMemory = FALSE;                               // If it's local data and the right size, we can keep memory.
125         }
126
127         dest->rows = orig->rows;
128         dest->cols = orig->cols;
129         dest->colMajor = orig->colMajor;
130         dest->originalRows = dest->rows;
131         dest->originalCols = dest->cols;
132         dest->originalColMajor = dest->colMajor;
133         if (copyState) dest->currentState = orig->currentState;
134         dest->lastCompute = orig->lastCompute;
135         dest->lastRow = orig->lastRow;
136         dest->hasMatrixNumber = orig->hasMatrixNumber;
137         dest->matrixNumber = orig->matrixNumber;
138
139         if(dest->rows == 0 || dest->cols == 0) {
140                 omxFreeMatrixData(dest);
141                 dest->data = NULL;
142                 dest->localData=FALSE;
143         } else {
144                 if(regenerateMemory) {
145                         omxFreeMatrixData(dest);                                                                                        // Free and regenerate memory
146                         dest->data = (double*) Calloc(dest->rows * dest->cols, double);
147                 }
148                 memcpy(dest->data, orig->data, dest->rows * dest->cols * sizeof(double));
149                 dest->localData = TRUE;
150         }
151
152         dest->aliasedPtr = NULL;
153
154         omxMatrixCompute(dest);
155
156 }
157
158 void omxAliasMatrix(omxMatrix *dest, omxMatrix *src) {
159         omxCopyMatrix(dest, src, TRUE);
160         dest->aliasedPtr = src;                                 // Alias now follows back matrix precisely.
161 }
162
163 void omxFreeMatrixData(omxMatrix * om) {
164
165         if(om->localData && om->data != NULL) {
166                 if(OMX_DEBUG_MATRIX) { Rprintf("Freeing matrix at 0x%0x. Localdata = %d.\n", om->data, om->localData); }
167                 Free(om->data);
168                 om->data = NULL;
169                 om->localData = FALSE;
170         }
171
172 }
173
174 void omxFreeAllMatrixData(omxMatrix *om) {
175     
176     if(om == NULL) return;
177
178         if(OMX_DEBUG) { 
179             Rprintf("Freeing matrix at 0x%0x with data = 0x%x, algebra 0x%x, and objective 0x%x.\n", 
180                 om, om->data, om->algebra); 
181         }
182
183         if(om->localData && om->data != NULL) {
184                 Free(om->data);
185                 om->data = NULL;
186                 om->localData = FALSE;
187         }
188
189         if(om->algebra != NULL) {
190                 omxFreeAlgebraArgs(om->algebra);
191                 om->algebra = NULL;
192         }
193
194         if(om->objective != NULL) {
195                 omxFreeObjectiveArgs(om->objective);
196                 om->objective = NULL;
197         }
198         
199         if(om->isTemporary) {
200                 Free(om);
201                 om = NULL;
202         }
203
204 }
205
206 void omxZeroByZeroMatrix(omxMatrix *om) {
207         if (om->rows > 0 || om->cols > 0) {
208                 omxResizeMatrix(om, 0, 0, FALSE);
209         }
210 }
211
212 omxMatrix* omxNewIdentityMatrix(int nrows, omxState* state) {
213         omxMatrix* newMat = NULL;
214         int l,k;
215
216         newMat = omxInitMatrix(newMat, nrows, nrows, FALSE, state);
217         for(k =0; k < newMat->rows; k++) {
218                 for(l = 0; l < newMat->cols; l++) {
219                         if(l == k) {
220                                 omxSetMatrixElement(newMat, k, l, 1);
221                         } else {
222                                 omxSetMatrixElement(newMat, k, l, 0);
223                         }
224                 }
225         }
226         return newMat;
227 }
228
229 omxMatrix* omxDuplicateMatrix(omxMatrix* tgt, omxMatrix* src, omxState* newState, short fullCopy) {
230     omxMatrix* newMat;
231     
232     if(src == NULL) return NULL;
233     newMat = omxInitMatrix(tgt, src->rows, src->cols, FALSE, newState);
234     omxCopyMatrix(newMat, src, FALSE);
235     if(src->algebra != NULL) {
236         omxDuplicateAlgebraMatrix(newMat, src, newState, fullCopy);
237     }
238     if(src->objective != NULL) {
239         omxDuplicateObjectiveMatrix(newMat, src, newState, fullCopy);
240     }
241     
242     return newMat;
243     
244 }
245
246 void omxResizeMatrix(omxMatrix *om, int nrows, int ncols, unsigned short keepMemory) {
247         // Always Recompute() before you Resize().
248         if(OMX_DEBUG_MATRIX) { 
249                 Rprintf("Resizing matrix from (%d, %d) to (%d, %d) (keepMemory: %d)", 
250                         om->rows, om->cols, 
251                         nrows, ncols, keepMemory);
252         }
253         if((keepMemory == FALSE) && (om->rows != nrows || om->cols != ncols)) {
254                 if(OMX_DEBUG_MATRIX) { Rprintf(" and regenerating memory to do it"); }
255                 omxFreeMatrixData(om);
256                 om->data = (double*) Calloc(nrows * ncols, double);
257                 om->localData = TRUE;
258         } else if(om->originalRows * om->originalCols < nrows * ncols) {
259                 warning("Upsizing an existing matrix may cause undefined behavior.\n"); // TODO: Define this behavior?
260         }
261
262         if(OMX_DEBUG_MATRIX) { Rprintf(".\n"); }
263         om->rows = nrows;
264         om->cols = ncols;
265         if(keepMemory == FALSE) {
266                 om->originalRows = om->rows;
267                 om->originalCols = om->cols;
268         }
269
270         omxMatrixCompute(om);
271 }
272
273 void omxResetAliasedMatrix(omxMatrix *om) {
274         om->rows = om->originalRows;
275         om->cols = om->originalCols;
276         if(om->aliasedPtr != NULL) {
277                 memcpy(om->data, om->aliasedPtr->data, om->rows*om->cols*sizeof(double));
278                 om->colMajor = om->aliasedPtr->colMajor;
279         }
280         omxMatrixCompute(om);
281 }
282
283 void omxMatrixCompute(omxMatrix *om) {
284
285         if(OMX_DEBUG_MATRIX) { Rprintf("Matrix compute: 0x%0x, 0x%0x, algebra: 0x%x.\n", om, om->currentState, om->algebra); }
286         om->majority = &(omxMatrixMajorityList[(om->colMajor?1:0)]);
287         om->minority = &(omxMatrixMajorityList[(om->colMajor?0:1)]);
288         om->leading = (om->colMajor?om->rows:om->cols);
289         om->lagging = (om->colMajor?om->cols:om->rows);
290
291         for(int i = 0; i < om->numPopulateLocations; i++) {
292                 omxRecompute(om->populateFrom[i]);                              // Make sure it's up to date
293                 double value = omxMatrixElement(om->populateFrom[i], om->populateFromRow[i], om->populateFromCol[i]);
294                 omxSetMatrixElement(om, om->populateToRow[i], om->populateToCol[i], value);
295                 // And then fill in the details.  Use the Helper here in case of transposition/downsampling.
296         }
297
298         om->isDirty = FALSE;
299         om->lastCompute = om->currentState->computeCount;
300         om->lastRow = om->currentState->currentRow;
301
302 }
303
304 double* omxLocationOfMatrixElement(omxMatrix *om, int row, int col) {
305         int index = 0;
306         if(om->colMajor) {
307                 index = col * om->rows + row;
308         } else {
309                 index = row * om->cols + col;
310         }
311         return om->data + index;
312 }
313
314 void vectorElementError(int index, int numrow, int numcol) {
315         char *errstr = calloc(250, sizeof(char));
316         sprintf(errstr, "Requested improper index (%d) from (%d, %d) vector.", 
317                 index, numrow, numcol);
318         error(errstr);
319         free(errstr);
320 }
321
322 void setMatrixError(int row, int col, int numrow, int numcol) {
323         char *errstr = calloc(250, sizeof(char));
324         sprintf(errstr, "Attempted to set improper value (%d, %d) into (%d, %d) matrix.", 
325                 row, col, numrow, numcol);
326         error(errstr);
327         free(errstr);
328 }
329
330 void matrixElementError(int row, int col, int numrow, int numcol) {
331         char *errstr = calloc(250, sizeof(char));
332         sprintf(errstr, "Requested improper value (%d, %d) from (%d, %d) matrix.",
333                 row, col, numrow, numcol);
334         error(errstr);
335         free(errstr);
336 }
337
338 void setVectorError(int index, int numrow, int numcol) {
339         char *errstr = calloc(250, sizeof(char));
340         sprintf(errstr, "Setting improper index (%d) from (%d, %d) vector.", 
341                 index, numrow, numcol);
342         error(errstr);
343         free(errstr);
344 }
345
346 double omxAliasedMatrixElement(omxMatrix *om, int row, int col) {
347         int index = 0;
348         if(row >= om->originalRows || col >= om->originalCols) {
349                 char *errstr = calloc(250, sizeof(char));
350                 sprintf(errstr, "Requested improper value (%d, %d) from (%d, %d) matrix.", row+1, col+1, om->originalRows, om->originalCols);
351                 error(errstr);
352                 free(errstr);
353         return (NA_REAL);
354         }
355         if(om->colMajor) {
356                 index = col * om->originalRows + row;
357         } else {
358                 index = row * om->originalCols + col;
359         }
360         return om->data[index];
361 }
362
363 void omxMarkDirty(omxMatrix *om) { om->isDirty = TRUE; }
364
365 unsigned short omxMatrixNeedsUpdate(omxMatrix *om) {
366         for(int i = 0; i < om->numPopulateLocations; i++) {
367                 if(omxNeedsUpdate(om->populateFrom[i])) return TRUE;    // Make sure it's up to date
368         }
369     return FALSE;
370 };
371
372 omxMatrix* omxNewMatrixFromMxMatrix(SEXP mxMatrix, omxState* state) {
373 /* Creates and populates an omxMatrix with details from an R MxMatrix object. */
374         omxMatrix *om = NULL;
375         om = omxInitMatrix(NULL, 0, 0, FALSE, state);
376         return omxFillMatrixFromMxMatrix(om, mxMatrix, state);
377 }
378
379 omxMatrix* omxNewMatrixFromRPrimitive(SEXP rObject, omxState* state, 
380         unsigned short hasMatrixNumber, int matrixNumber) {
381 /* Creates and populates an omxMatrix with details from an R matrix object. */
382         omxMatrix *om = NULL;
383         om = omxInitMatrix(NULL, 0, 0, FALSE, state);
384         return omxFillMatrixFromRPrimitive(om, rObject, state, hasMatrixNumber, matrixNumber);
385 }
386
387 omxMatrix* omxFillMatrixFromMxMatrix(omxMatrix* om, SEXP mxMatrix, omxState* state) {
388 /* Populates the fields of a omxMatrix with details from an R Matrix. */
389         if(inherits(mxMatrix, "MxMatrix")) {
390                 if(OMX_DEBUG) { Rprintf("R matrix is Mx Matrix.  Processing.\n"); }
391                 SEXP matrix;
392                 PROTECT(matrix = GET_SLOT(mxMatrix,  install("values")));
393                 om = fillMatrixHelperFunction(om, matrix, state, 0, 0);
394                 UNPROTECT(1);
395         } else {
396                 error("Recieved unknown matrix type in omxFillMatrixFromMxMatrix.");
397         }
398         return(om);
399 }
400
401 omxMatrix* omxFillMatrixFromRPrimitive(omxMatrix* om, SEXP rObject, omxState* state,
402         unsigned short hasMatrixNumber, int matrixNumber) {
403 /* Populates the fields of a omxMatrix with details from an R object. */
404         if(!isMatrix(rObject) && !isVector(rObject)) { // Sanity Check
405                 error("Recieved unknown matrix type in omxFillMatrixFromRPrimitive.");
406         }
407         return(fillMatrixHelperFunction(om, rObject, state, hasMatrixNumber, matrixNumber));
408 }
409
410
411
412 omxMatrix* fillMatrixHelperFunction(omxMatrix* om, SEXP matrix, omxState* state,
413         unsigned short hasMatrixNumber, int matrixNumber) {
414
415         SEXP matrixDims;
416         int* dimList;
417
418         if(OMX_DEBUG) { Rprintf("Filling omxMatrix from R matrix.\n"); }
419
420         if(om == NULL) {
421                 om = omxInitMatrix(NULL, 0, 0, FALSE, state);
422         }
423
424         om->data = REAL(AS_NUMERIC(matrix));    // TODO: Class-check first?
425
426         if(isMatrix(matrix)) {
427                 PROTECT(matrixDims = getAttrib(matrix, R_DimSymbol));
428                 dimList = INTEGER(matrixDims);
429                 om->rows = dimList[0];
430                 om->cols = dimList[1];
431                 UNPROTECT(1);   // MatrixDims
432         } else if (isVector(matrix)) {          // If it's a vector, assume it's a row vector. BLAS doesn't care.
433                 if(OMX_DEBUG) { Rprintf("Vector discovered.  Assuming rowity.\n"); }
434                 om->rows = 1;
435                 om->cols = length(matrix);
436         }
437         if(OMX_DEBUG) { Rprintf("Matrix connected to (%d, %d) matrix or MxMatrix.\n", om->rows, om->cols); }
438
439         om->localData = FALSE;
440         om->colMajor = TRUE;
441         om->originalRows = om->rows;
442         om->originalCols = om->cols;
443         om->originalColMajor = TRUE;
444         om->aliasedPtr = NULL;
445         om->algebra = NULL;
446         om->objective = NULL;
447         om->currentState = state;
448         om->lastCompute = -1;
449         om->lastRow = -1;
450         om->hasMatrixNumber = hasMatrixNumber;
451         om->matrixNumber = matrixNumber;
452
453
454         if(OMX_DEBUG) { Rprintf("Pre-compute call.\n");}
455         omxMatrixCompute(om);
456         if(OMX_DEBUG) { Rprintf("Post-compute call.\n");}
457
458         if(OMX_DEBUG) {
459                 omxPrintMatrix(om, "Finished importing matrix");
460         }
461
462         return om;
463 }
464
465 void omxProcessMatrixPopulationList(omxMatrix* matrix, SEXP matStruct) {
466
467         if(OMX_DEBUG) { Rprintf("Processing Population List: %d elements.\n", length(matStruct) - 1); }
468         SEXP subList;
469
470         if(length(matStruct) > 1) {
471                 int numPopLocs = length(matStruct) - 1;
472                 matrix->numPopulateLocations = numPopLocs;
473                 matrix->populateFrom = (omxMatrix**)R_alloc(numPopLocs, sizeof(omxMatrix*));
474                 matrix->populateFromRow = (int*)R_alloc(numPopLocs, sizeof(int));
475                 matrix->populateFromCol = (int*)R_alloc(numPopLocs, sizeof(int));
476                 matrix->populateToRow = (int*)R_alloc(numPopLocs, sizeof(int));
477                 matrix->populateToCol = (int*)R_alloc(numPopLocs, sizeof(int));
478         }
479
480         for(int i = 0; i < length(matStruct)-1; i++) {
481                 PROTECT(subList = AS_INTEGER(VECTOR_ELT(matStruct, i+1)));
482
483                 int* locations = INTEGER(subList);
484                 int loc = locations[0];
485                 if(OMX_DEBUG) { Rprintf("."); } //:::
486                 if(loc < 0) {                   // NOTE: This duplicates some of the functionality of NewMatrixFromMxIndex
487                         matrix->populateFrom[i] = matrix->currentState->matrixList[(~loc)];
488                 } else {
489                         matrix->populateFrom[i] = matrix->currentState->algebraList[(loc)];
490                 }
491                 matrix->populateFromRow[i] = locations[1];
492                 matrix->populateFromCol[i] = locations[2];
493                 matrix->populateToRow[i] = locations[3];
494                 matrix->populateToCol[i] = locations[4];
495
496                 UNPROTECT(1); // subList
497         }
498 }
499
500 void omxToggleRowColumnMajor(omxMatrix *mat) {
501
502         int i, j;
503         int nrows = mat->rows;
504         int ncols = mat->cols;
505         
506         double *newdata = (double*) Calloc(nrows * ncols, double);
507         double *olddata = mat->data;
508
509         if (mat->colMajor) {
510                 for(i = 0; i < ncols; i++) {
511                         for(j = 0; j < nrows; j++) {
512                                 newdata[i + ncols * j] = olddata[i * nrows + j];
513                         }
514                 }
515         } else {
516                 for(i = 0; i < nrows; i++) {
517                         for(j = 0; j < ncols; j++) {
518                                 newdata[i + nrows * j] = olddata[i * ncols + j];
519                         }
520                 }
521         }
522
523         if (mat->localData) {
524                 Free(mat->data);
525         }
526
527         mat->localData = TRUE;
528         mat->data = newdata;
529         mat->colMajor = !mat->colMajor;
530 }
531
532 void omxTransposeMatrix(omxMatrix *mat) {
533         mat->colMajor = !mat->colMajor;
534         omxMatrixCompute(mat);
535 }
536
537 void omxRemoveRowsAndColumns(omxMatrix *om, int numRowsRemoved, int numColsRemoved, int rowsRemoved[], int colsRemoved[])
538 {
539     // TODO: Create short-circuit form of omxRemoveRowsAndCols to remove just rows or just columns.
540 //      if(OMX_DEBUG_MATRIX) { Rprintf("Removing %d rows and %d columns from 0x%0x.\n", numRowsRemoved, numColsRemoved, om);}
541
542         if(numRowsRemoved < 1 && numColsRemoved < 1) { return; }
543
544         int oldRows, oldCols;
545
546         if(om->aliasedPtr == NULL) {
547                 if(om->originalRows == 0 || om->originalCols == 0) {
548                         om->originalRows = om->rows;
549                         om->originalCols = om->cols;
550                 }
551                 oldRows = om->originalRows;
552                 oldCols = om->originalCols;
553         } else {
554                 oldRows = om->aliasedPtr->rows;
555                 oldCols = om->aliasedPtr->cols;
556         }
557
558         int nextCol = 0;
559         int nextRow = 0;
560
561         if(om->rows > om->originalRows || om->cols > om->originalCols) {        // sanity check.
562                 error("Aliased Matrix is too small for alias.");
563         }
564
565         om->rows = oldRows - numRowsRemoved;
566         om->cols = oldCols - numColsRemoved;
567
568         // Note:  This really aught to be done using a matrix multiply.  Why isn't it?
569         for(int j = 0; j < oldCols; j++) {
570                 if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Handling column %d/%d...", j, oldCols);}
571                 if(colsRemoved[j]) {
572                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Removed.\n");}
573                         continue;
574                 } else {
575                         nextRow = 0;
576                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Rows (max %d): ", oldRows); }
577                         for(int k = 0; k < oldRows; k++) {
578                                 if(rowsRemoved[k]) {
579                                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("%d removed....", k);}
580                                         continue;
581                                 } else {
582                                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("%d kept....", k);}
583                                         if(om->aliasedPtr == NULL) {
584                                                 if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Self-aliased matrix access.\n");}
585                                                 omxSetMatrixElement(om, nextRow, nextCol, omxAliasedMatrixElement(om, k, j));
586                                         } else {
587                                                 if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("Matrix 0x%x re-aliasing to 0x%x.\n", om, om->aliasedPtr);}
588                                                 omxSetMatrixElement(om, nextRow, nextCol, omxMatrixElement(om->aliasedPtr, k,  j));
589                                         }
590                                         nextRow++;
591                                 }
592                         }
593                         if(OMX_DEBUG_MATRIX || OMX_DEBUG_ALGEBRA) { Rprintf("\n");}
594                         nextCol++;
595                 }
596         }
597
598         omxMatrixCompute(om);
599 }
600
601 /* Function wrappers that switch based on inclusion of algebras */
602 void omxPrint(omxMatrix *source, char* d) {                                     // Pretty-print a (small) matrix
603         if(source->algebra != NULL) omxAlgebraPrint(source->algebra, d);
604         else if(source->objective != NULL) omxObjectivePrint(source->objective, d);
605         else omxPrintMatrix(source, d);
606 }
607
608 unsigned short omxNeedsUpdate(omxMatrix *matrix) {
609         unsigned short retval;
610         /* Simplest update check: If we're dirty or haven't computed this cycle (iteration or row), we need to. */
611         // TODO : Implement a dependency-tree-based dirtiness propagation system
612         if(OMX_DEBUG_MATRIX) {Rprintf("Matrix 0x%x NeedsUpdate?", matrix);}
613
614         if(matrix == NULL) {
615                 if(OMX_DEBUG_MATRIX) {Rprintf("matrix argument is NULL. ");}
616                 retval = FALSE;         // Not existing means never having to say you need to recompute.
617         } else if(matrix->isDirty) {
618                 if(OMX_DEBUG_MATRIX) {Rprintf("matrix is dirty. ");}
619                 retval = TRUE;
620         } else if(matrix->lastCompute < matrix->currentState->computeCount) {
621                 if(OMX_DEBUG_MATRIX) {Rprintf("matrix last compute is less than current compute count. ");}
622                 retval = TRUE;          // No need to check args if oa's dirty.
623         } else if(matrix->lastRow != matrix->currentState->currentRow) {
624                 if(OMX_DEBUG_MATRIX) {Rprintf("matrix last row is less than current row. ");}
625                 retval = TRUE;                  // Ditto.
626         } else if(matrix->algebra != NULL) {
627                 if(OMX_DEBUG_MATRIX) {Rprintf("checking algebra needs update. ");}
628                 retval = omxAlgebraNeedsUpdate(matrix->algebra);
629         } else if(matrix->objective != NULL) {
630                 if(OMX_DEBUG_MATRIX) {Rprintf("checking objective function needs update. ");}
631                 retval = omxObjectiveNeedsUpdate(matrix->objective);
632         } else {
633                 if(OMX_DEBUG_MATRIX) {Rprintf("checking matrix needs update. ");}
634                 retval = omxMatrixNeedsUpdate(matrix);
635         }
636         if(OMX_DEBUG_MATRIX && retval) {Rprintf("Yes.\n");}
637         if(OMX_DEBUG_MATRIX && !retval) {Rprintf("No.\n");}
638         return(retval);
639 }
640
641 void omxRecompute(omxMatrix *matrix) {
642         if(!omxNeedsUpdate(matrix)) return;
643         if(matrix->algebra != NULL) omxAlgebraCompute(matrix->algebra);
644         else if(matrix->objective != NULL) omxObjectiveCompute(matrix->objective);
645         else omxMatrixCompute(matrix);
646 }
647
648 void omxCompute(omxMatrix *matrix) {
649         if(matrix->algebra != NULL) omxAlgebraCompute(matrix->algebra);
650         else if(matrix->objective != NULL) omxObjectiveCompute(matrix->objective);
651         else omxMatrixCompute(matrix);
652 }