EpetraExt_MatlabEngine.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //                     New_Package Example Package
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #include "EpetraExt_MatlabEngine.h"
00030 #include "EpetraExt_PutMultiVector.h"
00031 #include "EpetraExt_PutRowMatrix.h"
00032 #include "EpetraExt_PutBlockMap.h"
00033 
00034 #include "Epetra_Comm.h"
00035 #include "Epetra_MultiVector.h"
00036 #include "Epetra_SerialDenseMatrix.h"
00037 #include "Epetra_IntSerialDenseMatrix.h"
00038 #include "Epetra_IntVector.h"
00039 #include "Epetra_RowMatrix.h"
00040 #include "Epetra_DataAccess.h"
00041 #include "Epetra_Import.h"
00042 #include "Epetra_Export.h"
00043 #include "Epetra_CombineMode.h"
00044 #include "Epetra_CrsMatrix.h"
00045 #include "Epetra_Map.h"
00046 #include "Epetra_CrsMatrix.h"
00047 
00048 using namespace EpetraExt;
00049 namespace EpetraExt {
00050 
00051 //=============================================================================
00052 EpetraExt_MatlabEngine::EpetraExt_MatlabEngine (const Epetra_Comm& Comm):Comm_(Comm) {
00053   if (Comm_.MyPID() == 0) {
00054     // construct the MATLAB engine
00055     Engine_ = engOpen(NULL);
00056   }
00057 } 
00058  
00059 //=============================================================================
00060 EpetraExt_MatlabEngine::~EpetraExt_MatlabEngine (void) {
00061   if (Comm_.MyPID() == 0) {
00062     // to destruct the MATLAB engine
00063     engClose(Engine_);
00064   }
00065 }
00066 
00067 //=============================================================================
00068 int EpetraExt_MatlabEngine::EvalString (char* command, char* outputBuffer, int outputBufferSize) {
00069   // send a string command to the MATLAB engine
00070   if (Comm_.MyPID() == 0) {
00071     if (outputBuffer != NULL) {
00072       if (engOutputBuffer(Engine_, outputBuffer, outputBufferSize)) {
00073         return(-4);
00074       }
00075     }
00076     if (engEvalString(Engine_, command)) {
00077       return(-3);
00078     }
00079   }
00080   
00081   return(0);
00082 }
00083 
00084 //=============================================================================
00085 int EpetraExt_MatlabEngine::PutMultiVector(const Epetra_MultiVector& A, const char * variableName) {
00086   mxArray* matlabA = 0;
00087   double* pr = 0;
00088   if (Comm_.MyPID() == 0) {
00089     matlabA = mxCreateDoubleMatrix(A.GlobalLength(), A.NumVectors(), mxREAL);
00090     pr = mxGetPr(matlabA);
00091   }
00092 
00093   if (Matlab::CopyMultiVector(&pr, A)) {
00094     mxDestroyArray(matlabA);
00095     return(-2);
00096   }
00097 
00098   if (Comm_.MyPID() == 0)
00099     if (PutIntoMatlab(variableName, matlabA)) {
00100     mxDestroyArray(matlabA);
00101     return(-1);
00102   }
00103   
00104   mxDestroyArray(matlabA);
00105   return(0);
00106 }
00107 
00108 //=============================================================================
00109 int EpetraExt_MatlabEngine::PutRowMatrix(const Epetra_RowMatrix& A, const char* variableName, bool transA) {
00110   mxArray* matlabA = 0;
00111   if (Comm_.MyPID() == 0) {
00112     // since matlab uses column major for matrices, switch row and column numbers
00113     matlabA = mxCreateSparse(A.OperatorDomainMap().MaxAllGID() - A.OperatorDomainMap().MinAllGID()+1, A.OperatorRangeMap().MaxAllGID() - A.OperatorRangeMap().MinAllGID() + 1, A.NumGlobalNonzeros(), mxREAL);
00114   }
00115   //cout << "numrows: " << A.RowMatrixColMap().NumGlobalElements() << " numCols: " << A.NumGlobalRows() << "numNonZeros: " << A.NumGlobalNonzeros() << "\n";
00116 
00117   //cout << "calling CopyRowMatrix\n";
00118   if (Matlab::CopyRowMatrix(matlabA, A)) {
00119     mxDestroyArray(matlabA);
00120     return(-2);
00121   }
00122 
00123   //cout << "done doing CopyRowMatrix\n";
00124   if (Comm_.MyPID() == 0) {
00125 
00126     /*cout << "printing matlabA pointers\n";
00127     double* matlabAvaluesPtr = mxGetPr(matlabA);
00128     int* matlabAcolumnIndicesPtr = mxGetJc(matlabA);
00129     int* matlabArowIndicesPtr = mxGetIr(matlabA);
00130     for(int i=0; i < A.NumGlobalNonzeros(); i++) {
00131       cout << "*matlabAvaluesPtr: " << *matlabAvaluesPtr++ << " *matlabAcolumnIndicesPtr: " << *matlabAcolumnIndicesPtr++ << " *matlabArowIndicesPtr" << *matlabArowIndicesPtr++ << "\n";
00132     }
00133     cout << "done printing matlabA pointers\n";
00134     */
00135     if (PutIntoMatlab(variableName, matlabA)) {
00136     mxDestroyArray(matlabA);
00137     return(-1);
00138     }
00139 
00140     if (!transA) {
00141       char* buff = new char[128];;
00142       sprintf(buff, "%s = %s'", variableName, variableName);
00143       if (EvalString(buff)) {
00144         mxDestroyArray(matlabA);
00145         return(-3);
00146       }
00147     }
00148   }
00149 
00150   //cout << "done with everything in PutRowMatrix, going to destroy matlabA\n" << "matlabA=" << matlabA << "\n";
00151   mxDestroyArray(matlabA);
00152   //cout << "done destroying matlabA\n";
00153   return(0);
00154 }
00155 
00156 //=============================================================================
00157 int EpetraExt_MatlabEngine::PutCrsGraph(const Epetra_CrsGraph& A, const char* variableName, bool transA) {
00158   return(-9999);
00159 }
00160 
00161 //=============================================================================
00162 int EpetraExt_MatlabEngine::PutSerialDenseMatrix(const Epetra_SerialDenseMatrix& A, const char* variableName, int proc) {
00163   if (proc == 0) {
00164     if (Comm_.MyPID() == 0) {
00165       mxArray* matlabA = 0;
00166 
00167       int numRows = A.M();
00168       int numCols = A.N();
00169 
00170       matlabA = mxCreateDoubleMatrix(numRows, numCols, mxREAL);
00171 
00172       int row;
00173       int col;
00174       double* targetPtr = 0;
00175       double* sourcePtr = 0;
00176       double* source = (double*)A.A();
00177       double* target = (double*)mxGetPr(matlabA);
00178       int source_LDA = A.LDA();
00179       int target_LDA = A.LDA();
00180       for (col = 0; col < numCols; col++) {
00181         targetPtr = target + (col * target_LDA);
00182         sourcePtr = source + (col * source_LDA);
00183         for (row = 0; row < numRows; row++) {
00184           *targetPtr++ = *sourcePtr++;
00185         }
00186       }
00187 
00188       if (PutIntoMatlab(variableName, matlabA)) {
00189         mxDestroyArray(matlabA);
00190         return(-1);
00191       }
00192 
00193       mxDestroyArray(matlabA);
00194     }
00195 
00196     return(0);
00197   }
00198   else {
00199     Epetra_MultiVector* tempMultiVector = 0;
00200     if (proc == Comm_.MyPID()) {
00201       int* numVectors = new int[1];
00202       int* temp = new int[1];
00203       temp[0] = A.N();
00204       Comm_.MaxAll (temp, numVectors, 1);
00205       const Epetra_BlockMap tempBlockMap (-1, A.LDA(), 1, 0, Comm_);
00206       tempMultiVector = new Epetra_MultiVector (View, tempBlockMap, A.A(), A.LDA(), A.N());
00207     }
00208     else {
00209       int* numVectors = new int[1];
00210       int* temp = new int[1];
00211       temp[0] = 0;
00212       Comm_.MaxAll (temp, numVectors, 1);
00213       const Epetra_BlockMap tempBlockMap (-1, 0, 1, 0, Comm_);
00214       tempMultiVector = new Epetra_MultiVector (tempBlockMap, numVectors[0], false);
00215     }
00216 
00217     return(PutMultiVector(*tempMultiVector, variableName)); 
00218   }
00219 }
00220 
00221 //=============================================================================
00222 int EpetraExt_MatlabEngine::PutIntSerialDenseMatrix(const Epetra_IntSerialDenseMatrix& A, const char* variableName, int proc) {
00223   mxArray* matlabA = 0;
00224   
00225   if (proc == 0) {
00226     if (Comm_.MyPID() == 0) {
00227       int numRows = A.M();
00228       int numCols = A.N();
00229 
00230       matlabA = mxCreateDoubleMatrix(numRows, numCols, mxREAL);
00231 
00232       int row;
00233       int col;
00234       double* targetPtr = 0;
00235       int* sourcePtr = 0;
00236       int* source = (int*)A.A();
00237       double* target = (double*)mxGetPr(matlabA);
00238       int source_LDA = A.LDA();
00239       int target_LDA = A.LDA();
00240       for (col = 0; col < numCols; col++) {
00241         targetPtr = target + (col * target_LDA);
00242         sourcePtr = source + (col * source_LDA);
00243         for (row = 0; row < numRows; row++) {
00244           *targetPtr++ = *sourcePtr++;
00245         }
00246       }
00247 
00248       if (PutIntoMatlab(variableName, matlabA)) {
00249         mxDestroyArray(matlabA);
00250         return(-1);
00251       }
00252     }
00253   }
00254   else {
00255     int* numVectors = new int[2];
00256     if (proc == Comm_.MyPID()) {
00257       int* temp = new int[2];
00258       temp[0] = A.N();
00259       temp[1] = A.M();
00260       Comm_.MaxAll (temp, numVectors, 2);
00261     }
00262     else {
00263       int* temp = new int[2];
00264       temp[0] = 0;
00265       temp[1] = 0;
00266       Comm_.MaxAll (temp, numVectors, 2);
00267     }
00268 
00269     int numRows = numVectors[1];
00270     int numCols = numVectors[0];
00271     double* targetPtr = 0;
00272     int* sourcePtr = 0;
00273     int row;
00274     double* target = 0;
00275     const Epetra_BlockMap* tgBlockMap = 0;
00276     if (Comm_.MyPID() == 0) {
00277       matlabA = mxCreateDoubleMatrix(numRows, numCols, mxREAL);
00278       target = (double*)mxGetPr(matlabA);
00279       tgBlockMap = new Epetra_BlockMap(-1, numRows, 1, 0, Comm_);
00280     }
00281     else {
00282       tgBlockMap = new Epetra_BlockMap(-1, 0, 1, 0, Comm_);
00283     }
00284 
00285     const Epetra_BlockMap* srcBlockMap = 0;
00286     Epetra_IntVector* srcIntVector = 0;
00287     Epetra_IntVector tgIntVector (*tgBlockMap, false);
00288     if (proc == Comm_.MyPID()) {
00289       srcBlockMap = new Epetra_BlockMap(-1, A.LDA(), 1, 0, Comm_);
00290     }
00291     else {
00292       srcBlockMap = new Epetra_BlockMap(-1, 0, 1, 0, Comm_);
00293     }
00294 
00295     Epetra_Import importer (*tgBlockMap, *srcBlockMap);
00296 
00297     for(int i=0; i < numVectors[0]; i++) {
00298       if (proc == Comm_.MyPID()) {
00299         srcIntVector = new Epetra_IntVector(View, *srcBlockMap, (int*)A[i]);
00300       }
00301       else {
00302         srcIntVector = new Epetra_IntVector(*srcBlockMap, false);
00303       }
00304       // need to add some error checking for this!
00305       tgIntVector.Import(*srcIntVector, importer, Insert);
00306       if (Comm_.MyPID() == 0) {
00307         targetPtr = target + (i * numRows);
00308         sourcePtr = (int*)tgIntVector.Values();
00309         for (row = 0; row < numRows; row++) {
00310           *targetPtr++ = *sourcePtr++;
00311         }
00312       }
00313     }
00314 
00315     if (PutIntoMatlab(variableName, matlabA)) {
00316       mxDestroyArray(matlabA);
00317       return(-1);
00318     }
00319   }
00320   
00321   mxDestroyArray(matlabA);
00322   return(0);
00323 }
00324 
00325 //=============================================================================
00326 int EpetraExt_MatlabEngine::PutBlockMap(const Epetra_BlockMap& blockMap, const char* variableName, bool transA) {
00327   mxArray* matlabA = 0;
00328   if (Comm_.MyPID() == 0) {    
00329     int M = blockMap.NumGlobalElements();
00330     int N = 1;
00331   
00332     if (blockMap.MaxElementSize()>1) N = 2; // Non-trivial block map, store element sizes in second column
00333   
00334     matlabA = mxCreateSparse(N, M, M*N, mxREAL);
00335     int* matlabAcolumnIndicesPtr = mxGetJc(matlabA);
00336     matlabAcolumnIndicesPtr += M;
00337     *matlabAcolumnIndicesPtr = M*N; // set max cap.
00338   }
00339   
00340   if (Matlab::CopyBlockMap(matlabA, blockMap)) {
00341     mxDestroyArray(matlabA);
00342     return(-2);
00343   }
00344 
00345   if (Comm_.MyPID() == 0) {
00346     if (PutIntoMatlab(variableName, matlabA)) {
00347       mxDestroyArray(matlabA);
00348       return(-1);
00349     }
00350 
00351     if (!transA) {
00352       char* buff = new char[128];;
00353       sprintf(buff, "%s = %s'", variableName, variableName);
00354       if (EvalString(buff)) {
00355         mxDestroyArray(matlabA);
00356         return(-3);
00357       }
00358     }
00359   }
00360 
00361   mxDestroyArray(matlabA);
00362   return(0);
00363 }
00364 
00365 //=============================================================================
00366 int EpetraExt_MatlabEngine::PutIntoMatlab(const char* variableName, mxArray* matlabA) {
00367   if (Comm_.MyPID() != 0) {
00368     return(0);
00369   }
00370 #ifdef USE_ENGPUTARRAY
00371   // for matlab versions < 6.5 (release 13)
00372   mxSetName(matlabA, variableName);
00373   return(engPutArray(Engine_, matlabA));
00374 #else
00375   // for matlab versions >= 6.5 (release 13)
00376   return(engPutVariable(Engine_, variableName, matlabA));
00377 #endif
00378 }
00379 
00380 //=============================================================================
00381 int EpetraExt_MatlabEngine::GetMultiVector(const char* variableName, Epetra_MultiVector& A) {
00382   mxArray* matlabA = 0;
00383   int ierr = 0;
00384   if (ierr = GetmxArray(variableName, &matlabA)) {
00385     mxDestroyArray(matlabA);
00386     return(ierr);
00387   }
00388 
00389   bool isSparse = false;
00390   int numRows = 0;
00391   int numCols = 0;
00392   int numNonZeros = 0;
00393 
00394   if (ierr = GetmxArrayDimensions(matlabA, isSparse, numRows, numCols, numNonZeros)) {
00395     mxDestroyArray(matlabA);
00396     return(ierr);
00397   }
00398 
00399   if ((Comm_.MyPID() == 0) && isSparse) {
00400     mxDestroyArray(matlabA);
00401     return(-7);
00402   }
00403   
00404   // tempMap is nontrivial only on PE0
00405   Epetra_Map tempMap (-1, numRows, 0, Comm_);
00406   
00407   int numVectors = 0;
00408   Comm_.MaxAll (&numCols, &numVectors, 1);
00409   double* tempValues = 0;
00410   if (Comm_.MyPID() == 0) {
00411     tempValues = mxGetPr(matlabA);
00412   }
00413   Epetra_MultiVector tempMV (View, tempMap, tempValues, numRows, numVectors);
00414   Epetra_Export exporter (tempMap, A.Map());
00415   A.Export(tempMV, exporter, Insert);
00416   
00417   mxDestroyArray(matlabA);
00418   return(0);
00419 }
00420 
00421 //=============================================================================
00422 int EpetraExt_MatlabEngine::GetSerialDenseMatrix(const char* variableName, Epetra_SerialDenseMatrix& A, int proc) {
00423   int ierr = 0;
00424   int numMyGIDs = 0;
00425   if (Comm_.MyPID() == proc) {
00426     numMyGIDs = A.M();
00427   }
00428   Epetra_Map tempMap (-1, numMyGIDs, 0, Comm_);
00429   Epetra_MultiVector tempMV (View, tempMap, A.A(), numMyGIDs, A.N());
00430   
00431   if (ierr = GetMultiVector(variableName, tempMV)) {
00432     return(ierr);
00433   }
00434   
00435   if (Comm_.MyPID() == proc) {
00436     double* aValues = A.A();
00437     double* mvValues = tempMV.Values();
00438     for(int i=0; i < A.M() * A.N(); i++) {
00439       *aValues++ = *mvValues++;
00440     }
00441   }
00442   
00443   return(0);
00444 }
00445 
00446 //=============================================================================
00447 int EpetraExt_MatlabEngine::GetIntSerialDenseMatrix(const char* variableName, Epetra_IntSerialDenseMatrix& A, int proc) {
00448   int ierr = 0;
00449   int numMyGIDs = 0;
00450   double* values = 0;
00451   if (Comm_.MyPID() == proc) {
00452     numMyGIDs = A.M();
00453     int* aValues = A.A();
00454     values = new double[A.M() * A.N()];
00455     for (int i=0; i < A.M() * A.N(); i++) {
00456       *values++ = *aValues++;
00457     }
00458   }
00459   Epetra_Map tempMap (-1, numMyGIDs, 0, Comm_);
00460   Epetra_MultiVector tempMV (View, tempMap, values, numMyGIDs, A.N());
00461   
00462   if (ierr = GetMultiVector(variableName, tempMV)) {
00463     return(ierr);
00464   }
00465   
00466   if (Comm_.MyPID() == proc) {
00467     int* aValues = A.A();
00468     double* mvValues = tempMV.Values();
00469     for(int i=0; i < A.M() * A.N(); i++) {
00470       *aValues++ = *mvValues++;
00471     }
00472   }
00473   
00474   return(0);
00475 }
00476 
00477 //=============================================================================
00478 int EpetraExt_MatlabEngine::GetCrsMatrix(const char* inVariableName, Epetra_CrsMatrix& A, bool getTrans) {
00479   const char* variableName;
00480   
00481   if (!getTrans) {
00482     int inVariableNameLength = strlen(inVariableName);
00483     char* buff = new char[inVariableNameLength*2 + 11];
00484     sprintf(buff, "TRANS_%s = %s'", inVariableName, inVariableName);
00485     if (EvalString(buff)) {
00486       return(-3);
00487     }
00488     char* tempStr = new char[inVariableNameLength + 7];
00489     sprintf(tempStr, "TRANS_%s", inVariableName);
00490     variableName = tempStr;
00491   }
00492   else {
00493     variableName = inVariableName;
00494   }
00495   
00496   mxArray* matlabA = 0;
00497   int ierr = 0;
00498   if (ierr = GetmxArray(variableName, &matlabA)) {
00499     mxDestroyArray(matlabA);
00500     return(ierr);
00501   }
00502   
00503   if (!getTrans) {
00504     char* buff = new char[strlen(variableName) + 7];
00505     sprintf(buff, "clear %s", variableName);
00506     if (EvalString(buff)) {
00507       return(-3);
00508     }
00509   }
00510   
00511   bool isSparse = false;
00512   int numRows = 0;
00513   int numCols = 0;
00514   int numNonZeros = 0;
00515 
00516   // note that numCols and numRows are transposed on purpose here
00517   // we will treat the column major mxArray as a row major mxArray
00518   if (ierr = GetmxArrayDimensions(matlabA, isSparse, numCols, numRows, numNonZeros)) {
00519     mxDestroyArray(matlabA);
00520     return(ierr);
00521   }
00522 
00523   if ((Comm_.MyPID() == 0) && !isSparse) {
00524     mxDestroyArray(matlabA);
00525     return(-8);
00526   }
00527   
00528   // tempMap is nontrivial only on PE0
00529   Epetra_Map tempMap (-1, numRows, 0, Comm_);
00530   
00531   int numVectors = 0;
00532   Comm_.MaxAll (&numCols, &numVectors, 1);
00533   Epetra_CrsMatrix tempCRSM (View, tempMap, numVectors);
00534   if (Comm_.MyPID() == 0) {
00535     int* rowIndex = mxGetJc(matlabA);
00536     int* colIndex = mxGetIr(matlabA);
00537     double* values = mxGetPr(matlabA);
00538     int numCols = 0;
00539     for(int row = 0; row < numRows; row++) {
00540       numCols = *(rowIndex+1) - *rowIndex;
00541       if (numCols > 0) {
00542         int* indices = new int[numCols];
00543         int* indicesPtr = indices;
00544         for(int col=0; col < numCols; col++) {
00545           *indicesPtr++ = *colIndex++;
00546         }
00547         tempCRSM.InsertGlobalValues(row, numCols, values, indices);
00548       }
00549       values += numCols;
00550       rowIndex++;
00551     }
00552   }
00553   
00554   tempCRSM.FillComplete();
00555   Epetra_Export exporter (tempMap, A.RowMatrixRowMap());
00556   A.Export(tempCRSM, exporter, Insert);
00557   
00558   mxDestroyArray(matlabA);
00559   return(0);
00560 }
00561 
00562 //=============================================================================
00563 int EpetraExt_MatlabEngine::GetmxArrayDimensions(mxArray* matlabA, bool& isSparse, int& numRows, int& numCols, int& numNonZeros) {
00564   if (Comm_.MyPID() == 0) {
00565     if (mxGetNumberOfDimensions(matlabA) > 2) {
00566       return(-6);
00567     }
00568 
00569     const int* dimensions = mxGetDimensions(matlabA);
00570     numRows = dimensions[0];
00571     numCols = dimensions[1];
00572     isSparse = mxIsSparse(matlabA);
00573 
00574     if (isSparse) {
00575       numNonZeros = mxGetNzmax(matlabA);
00576     }
00577     else {
00578       numNonZeros = numRows * numCols;
00579     }
00580   }
00581   
00582   return(0);
00583 }
00584 
00585 //=============================================================================
00586 int EpetraExt_MatlabEngine::GetmxArray(const char* variableName, mxArray** matlabA) {
00587   if (Comm_.MyPID() == 0) {
00588 #ifdef USE_ENGPUTARRAY
00589   // for matlab versions < 6.5 (release 13)
00590   *matlabA = engGetArray(Engine_, variableName);
00591 #else
00592   // for matlab versions >= 6.5 (release 13)
00593   *matlabA = engGetVariable(Engine_, variableName);
00594 #endif
00595     if (matlabA == NULL) {
00596       return(-5);
00597     }
00598   }
00599   
00600   return(0);
00601 }
00602 
00603 } // namespace EpetraExt

Generated on Thu Sep 18 12:31:43 2008 for EpetraExt by doxygen 1.3.9.1