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