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