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