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