Epetra Package Browser (Single Doxygen Collection) Development
Epetra_OskiMatrix.cpp
Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2011 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 
00043 #include "Epetra_ConfigDefs.h"
00044 #include "Epetra_Map.h"
00045 
00046 #ifdef HAVE_OSKI
00047 #ifdef HAVE_EPETRA_TEUCHOS
00048 #include "Epetra_OskiMatrix.h"
00049 #include "Epetra_Import.h"
00050 
00051 //=============================================================================
00052 
00053 Epetra_OskiMatrix::Epetra_OskiMatrix(const Epetra_OskiMatrix& Source) 
00054   : Epetra_CrsMatrix(Source), 
00055   Epetra_View_(Source.Epetra_View_),
00056   Copy_Created_(true) {
00057     A_tunable_ = oski_CopyMat(Source.A_tunable_);
00058 }
00059 
00060 Epetra_OskiMatrix::Epetra_OskiMatrix(const Epetra_CrsMatrix& Source, 
00061              const Teuchos::ParameterList& List) 
00062   : Epetra_CrsMatrix(Source), 
00063   Epetra_View_(&Source) {
00064     bool AutoTune = false;
00065     bool DeepCopy = false;
00066     string Matrix = "general\0";
00067     bool IsDiagNotStored = false;
00068     bool IsArrayZeroBased = false;
00069     int MyIndexBase = 1;
00070     bool AreIndicesSorted = false;
00071     bool AreIndicesRepeated = false;
00072     oski_inmatprop_t MatrixType = MAT_GENERAL;
00073     oski_inmatprop_t Diagonal = MAT_DIAG_EXPLICIT;
00074     oski_inmatprop_t ArrayBasis = INDEX_ONE_BASED;
00075     oski_inmatprop_t SortedIndices = INDEX_UNSORTED;
00076     oski_inmatprop_t RepeatedIndices = INDEX_REPEATED;
00077     int* RowPtr = NULL;
00078     int* IndPtr = NULL;
00079     double* ValPtr = NULL;
00080     AutoTune = const_cast <Teuchos::ParameterList &>(List).get("autotune", false);
00081     if(List.isParameter("deepcopy")) 
00082       DeepCopy = Teuchos::getParameter<bool>(List, "deepcopy");
00083     if(AutoTune){  //Use parameters from the Epetra matrix to set as many fields as possible
00084       if(LowerTriangular())
00085         MatrixType = MAT_TRI_LOWER;
00086       if(UpperTriangular())
00087         MatrixType = MAT_TRI_UPPER;
00088       if(Sorted())
00089         SortedIndices = INDEX_SORTED;
00090       MyIndexBase = IndexBase();
00091       if(MyIndexBase == 0)
00092         ArrayBasis = INDEX_ZERO_BASED;
00093       else if(MyIndexBase == 1)
00094         ArrayBasis = INDEX_ONE_BASED;
00095       else if(!List.isParameter("zerobased")) {
00096         std::cerr << "An OskiMatrix must be either zero or one based.\n";
00097         return;
00098       }
00099       if(NoRedundancies())
00100         RepeatedIndices = INDEX_UNIQUE;
00101     }
00102     if(List.isParameter("matrixtype")) {
00103       Matrix = Teuchos::getParameter<string>(List, "matrixtype");
00104       if(Matrix == "general")
00105         MatrixType = MAT_GENERAL;
00106       else if(Matrix == "uppertri")
00107         MatrixType = MAT_TRI_UPPER;
00108       else if(Matrix == "lowertri")
00109         MatrixType = MAT_TRI_LOWER;
00110       else if(Matrix == "uppersymm")
00111         MatrixType = MAT_SYMM_UPPER;
00112       else if(Matrix == "lowersymm")
00113         MatrixType = MAT_SYMM_LOWER;
00114       else if(Matrix == "fullsymm")
00115         MatrixType = MAT_SYMM_FULL;
00116       else if(Matrix == "upperherm")
00117         MatrixType = MAT_HERM_UPPER;
00118       else if(Matrix == "lowerherm")
00119         MatrixType = MAT_HERM_LOWER;
00120       else if(Matrix == "fullherm")
00121         MatrixType = MAT_HERM_FULL;
00122     }
00123     if(List.isParameter("diagstored")) 
00124       IsDiagNotStored = Teuchos::getParameter<bool>(List, "diagstored");
00125     if(List.isParameter("zerobased")) 
00126       IsArrayZeroBased = Teuchos::getParameter<bool>(List, "zerobased");
00127     if(List.isParameter("sorted")) 
00128       AreIndicesSorted = Teuchos::getParameter<bool>(List, "sorted");
00129     if(List.isParameter("unique")) 
00130       DeepCopy = Teuchos::getParameter<bool>(List, "unique");
00131     if(IsDiagNotStored)
00132       Diagonal = MAT_UNIT_DIAG_IMPLICIT;
00133     if(IsArrayZeroBased)
00134       ArrayBasis = INDEX_ZERO_BASED;
00135     if(AreIndicesSorted)
00136       SortedIndices = INDEX_SORTED;
00137     if(AreIndicesRepeated)
00138       RepeatedIndices = INDEX_UNIQUE;
00139     if(ExtractCrsDataPointers(RowPtr, IndPtr, ValPtr)) {
00140       std::cerr << "Cannot wrap matrix as an Oski Matrix because at least one of FillComplete and Optimize Storage has not been called\n";
00141       return;
00142     }
00143     if(DeepCopy) {  
00144       Copy_Created_ = true; 
00145       A_tunable_ = oski_CreateMatCSR(RowPtr, IndPtr, ValPtr, Source.NumMyRows(), Source.NumMyCols(), COPY_INPUTMAT, 5, MatrixType, Diagonal, ArrayBasis, SortedIndices, RepeatedIndices);
00146     }
00147     else {
00148       Copy_Created_ = false;
00149       A_tunable_ = oski_CreateMatCSR(RowPtr, IndPtr, ValPtr, Source.NumMyRows(), Source.NumMyCols(), SHARE_INPUTMAT, 5, MatrixType, Diagonal, ArrayBasis, SortedIndices, RepeatedIndices);
00150     }
00151 }
00152 
00153 Epetra_OskiMatrix::~Epetra_OskiMatrix() {
00154   if(oski_DestroyMat(A_tunable_))
00155     std::cerr << "Destroy Matrix failed.\n";
00156 }
00157 
00158 int Epetra_OskiMatrix::ReplaceMyValues(int MyRow, 
00159                int NumEntries, 
00160                double* Values, 
00161                int* Indices) {
00162   int ReturnVal = 0;
00163   if (Copy_Created_) {
00164     for(int i = 0; i < NumEntries; i++) {
00165       ReturnVal = oski_SetMatEntry(A_tunable_, MyRow, Indices[i], Values[i]);
00166       if(ReturnVal)
00167   break;
00168     }
00169     if(!ReturnVal)
00170       ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->ReplaceMyValues(MyRow, NumEntries, Values, Indices);
00171   }
00172   else
00173     ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->ReplaceMyValues(MyRow, NumEntries, Values, Indices);
00174   if(ReturnVal)
00175     std::cerr << "Error in ReplaceMyValues\n";
00176   return ReturnVal;
00177 }
00178 
00179 int Epetra_OskiMatrix::SumIntoMyValues(int MyRow, 
00180                int NumEntries, 
00181                double* Values, 
00182                int* Indices) {
00183   int ReturnVal = 0;
00184   if (Copy_Created_) {
00185     for(int i = 0; i < NumEntries; i++) {
00186       ReturnVal = oski_SetMatEntry(A_tunable_, MyRow, Indices[i], oski_GetMatEntry(A_tunable_, MyRow, Indices[i]));
00187       if(ReturnVal)
00188   break;
00189     }
00190     if(!ReturnVal)
00191       ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->SumIntoMyValues(MyRow, NumEntries, Values, Indices);
00192   }
00193   else
00194     ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->SumIntoMyValues(MyRow, NumEntries, Values, Indices);
00195   if(ReturnVal)
00196     std::cerr << "Error in SumIntoMyValues\n";
00197   return ReturnVal;
00198 }
00199 
00200 int Epetra_OskiMatrix::ExtractDiagonalCopy(Epetra_Vector& Diagonal) const {
00201   int ReturnValue = 0;
00202   ReturnValue = Epetra_View_->ExtractDiagonalCopy(Diagonal);
00203   if (ReturnValue)
00204     std::cerr << "Error in ExtractDiagonalCopy\n";
00205   return ReturnValue;  
00206 }
00207 
00208 int Epetra_OskiMatrix::ReplaceDiagonalValues(const Epetra_OskiVector& Diagonal) {
00209   int ReturnVal = 0;
00210   if (Copy_Created_) {
00211     ReturnVal = oski_SetMatDiagValues(A_tunable_, 0, Diagonal.Oski_View());
00212     if(!ReturnVal)
00213       ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->ReplaceDiagonalValues(*Diagonal.Epetra_View());
00214   }
00215   else
00216     ReturnVal = const_cast <Epetra_CrsMatrix*> (Epetra_View_)->ReplaceDiagonalValues(*Diagonal.Epetra_View());
00217   if(ReturnVal)
00218     std::cerr << "Error in ReplaceDiagonalValues\n";
00219   return ReturnVal;
00220 }
00221 
00222 int Epetra_OskiMatrix::Multiply(bool TransA, 
00223         const Epetra_Vector& x, 
00224         Epetra_Vector& y) const {
00225   int ReturnVal;
00226   ReturnVal = this->Multiply(TransA, x, y, 1.0, 0.0);
00227   return ReturnVal;
00228 }
00229 
00230 int Epetra_OskiMatrix::Multiply(bool TransA, 
00231         const Epetra_Vector& x, 
00232         Epetra_Vector& y, 
00233         double Alpha, 
00234         double Beta) const {
00235   int ReturnVal;
00236  
00237   if(!Filled())
00238     EPETRA_CHK_ERR(-1); // Matrix must be filled.
00239 
00240   double* xp = (double*) x.Values();
00241   double* yp = (double*) y.Values();
00242 
00243   Epetra_Vector * xcopy = 0;
00244   if (&x==&y && Importer()==0 && Exporter()==0) {
00245     xcopy = new Epetra_Vector(x);
00246     xp = (double *) xcopy->Values();
00247   }
00248   UpdateImportVector(1); // Refresh import and output vectors if needed
00249   UpdateExportVector(1);
00250 
00251   if(!TransA) {
00252 
00253     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
00254     if(Importer() != 0) {
00255       EPETRA_CHK_ERR(ImportVector_->Import(x, *Importer(), Insert));
00256       xp = (double*) ImportVector_->Values();
00257     }
00258 
00259     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00260     if(Exporter() != 0) 
00261       yp = (double*) ExportVector_->Values();
00262     
00263 
00264     oski_vecview_t oskiX;
00265     oski_vecview_t oskiY;
00266     if(Importer() != 0) 
00267       oskiX = oski_CreateVecView(xp,ImportVector_->MyLength(),1);
00268     else               
00269       oskiX = oski_CreateVecView(xp,x.MyLength(),1);
00270     if(Exporter() != 0) 
00271       oskiY = oski_CreateVecView(yp,ExportVector_->MyLength(),1);
00272     else               
00273       oskiY = oski_CreateVecView(yp,y.MyLength(),1);
00274 
00275     //Do actual computation
00276     ReturnVal = oski_MatMult(A_tunable_, OP_NORMAL, Alpha, oskiX, Beta, oskiY);
00277 
00278     if(Exporter() != 0) {
00279       y.PutScalar(0.0); // Make sure target is zero
00280       EPETRA_CHK_ERR(y.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
00281     }
00282     // Handle case of rangemap being a local replicated map
00283     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
00284   } //if(!TransA)
00285   else { // Transpose operation
00286 
00287     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
00288     if(Exporter() != 0) {
00289       EPETRA_CHK_ERR(ExportVector_->Import(x, *Exporter(), Insert));
00290       xp = (double*) ExportVector_->Values();
00291     }
00292 
00293     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
00294     if(Importer() != 0) 
00295       yp = (double*) ImportVector_->Values();
00296 
00297     oski_vecview_t oskiX;
00298     oski_vecview_t oskiY;
00299     if(Importer() != 0) 
00300       oskiY = oski_CreateVecView(yp,ImportVector_->MyLength(),1);
00301     else               
00302       oskiY = oski_CreateVecView(yp,y.MyLength(),1);
00303     if(Exporter() != 0) 
00304       oskiX = oski_CreateVecView(xp,ExportVector_->MyLength(),1);
00305     else               
00306       oskiX = oski_CreateVecView(xp,x.MyLength(),1);
00307     
00308     // Do actual computation
00309     ReturnVal = oski_MatMult(A_tunable_, OP_TRANS, Alpha, oskiX, Beta, oskiY);
00310 
00311     if(Importer() != 0) {
00312       y.PutScalar(0.0); // Make sure target is zero
00313       EPETRA_CHK_ERR(y.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
00314     }
00315     // Handle case of rangemap being a local replicated map
00316     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
00317   }
00318   if(ReturnVal)
00319     std::cerr << "OskiVector multiply error\n";
00320   UpdateFlops(2 * NumGlobalNonzeros64());
00321   return ReturnVal;
00322 }
00323 
00324 int Epetra_OskiMatrix::Multiply(bool TransA, 
00325         const Epetra_MultiVector& X, 
00326         Epetra_MultiVector& Y) const {
00327   int ReturnVal;
00328   ReturnVal = this->Multiply(TransA, X, Y, 1.0, 0.0);
00329   return ReturnVal;
00330 }
00331 
00332 int Epetra_OskiMatrix::Multiply(bool TransA, 
00333         const Epetra_MultiVector& X, 
00334         Epetra_MultiVector& Y, 
00335         double Alpha, 
00336         double Beta) const {
00337   int ReturnVal;
00338 
00339   if(!Filled())
00340     EPETRA_CHK_ERR(-1); // Matrix must be filled.
00341 
00342   int NumVectors = X.NumVectors();
00343   if (NumVectors!=Y.NumVectors()) {
00344     EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
00345   }
00346 
00347   double** Xp = (double**) X.Pointers();
00348   double** Yp = (double**) Y.Pointers();
00349 
00350   int LDX = X.ConstantStride() ? X.Stride() : 0;
00351   int LDY = Y.ConstantStride() ? Y.Stride() : 0;
00352 
00353   Epetra_MultiVector* Xcopy = 0;
00354   if (&X==&Y && Importer()==0 && Exporter()==0) {
00355     Xcopy = new Epetra_MultiVector(X);
00356     Xp = (double **) Xcopy->Pointers();
00357     LDX = Xcopy->ConstantStride() ? Xcopy->Stride() : 0;
00358   }
00359   UpdateImportVector(NumVectors); // Refresh import and output vectors if needed
00360   UpdateExportVector(NumVectors);
00361 
00362   if (!TransA) {
00363 
00364     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
00365     if (Importer()!=0) {
00366       EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
00367       Xp = (double**)ImportVector_->Pointers();
00368       LDX = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
00369     }
00370 
00371     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00372     if (Exporter()!=0) {
00373       Yp = (double**)ExportVector_->Pointers();
00374       LDY = ExportVector_->ConstantStride() ? ExportVector_->Stride() : 0;
00375     }
00376 
00377     oski_vecview_t oskiX;
00378     oski_vecview_t oskiY;
00379     if (Importer()!=0) 
00380       oskiX = oski_CreateMultiVecView(*Xp,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
00381     else               
00382       oskiX = oski_CreateMultiVecView(*Xp,X.MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
00383     if (Exporter()!=0) 
00384       oskiY = oski_CreateMultiVecView(*Yp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
00385     else               
00386       oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
00387     // Do actual computation
00388     ReturnVal = oski_MatMult(A_tunable_, OP_NORMAL, Alpha, oskiX, Beta, oskiY);
00389     if(ReturnVal)
00390       std::cerr << "OskiMultiVector multiply error\n";
00391     if (Exporter()!=0) {
00392       Y.PutScalar(0.0);  // Make sure target is zero
00393       Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
00394     }
00395     // Handle case of rangemap being a local replicated map
00396     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
00397   }
00398   else { // Transpose operation
00399 
00400     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
00401 
00402     if (Exporter()!=0) {
00403       EPETRA_CHK_ERR(ExportVector_->Import(X, *Exporter(), Insert));
00404       Xp = (double**)ExportVector_->Pointers();
00405       LDX = ExportVector_->ConstantStride() ? ExportVector_->Stride() : 0;
00406     }
00407 
00408     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
00409     if (Importer()!=0) {
00410       Yp = (double**)ImportVector_->Pointers();
00411       LDY = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
00412     }
00413     
00414     oski_vecview_t oskiX;
00415     oski_vecview_t oskiY;
00416     if (Importer()!=0) 
00417       oskiY = oski_CreateMultiVecView(*Yp,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
00418     else               
00419       oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
00420     if (Exporter()!=0) 
00421       oskiX = oski_CreateMultiVecView(*Xp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
00422     else               
00423       oskiX = oski_CreateMultiVecView(*Xp,X.MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
00424 
00425     // Do actual computation
00426     ReturnVal = oski_MatMult(A_tunable_, OP_TRANS, Alpha, oskiX, Beta, oskiY);
00427     if(ReturnVal)
00428       std::cerr << "OskiMultiVector multiply error\n";
00429     if (Importer()!=0) {
00430       Y.PutScalar(0.0);  // Make sure target is zero
00431       EPETRA_CHK_ERR(Y.Export(*ImportVector_, *Importer(), Add)); // Fill Y with Values from export vector
00432     }
00433     // Handle case of rangemap being a local replicated map
00434     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1)  EPETRA_CHK_ERR(Y.Reduce());
00435   }
00436   UpdateFlops(2 * NumGlobalNonzeros64());
00437   //Y.ResetView(Yp);
00438   return ReturnVal;
00439 }
00440 
00441 int Epetra_OskiMatrix::Solve(bool Upper, bool TransA, bool UnitDiagonal, const Epetra_Vector& x, Epetra_Vector& y) const {
00442   int ReturnVal;
00443   ReturnVal = this->Solve(TransA, x, y, 1.0);
00444   return ReturnVal;
00445 }
00446 
00447 int Epetra_OskiMatrix::Solve(bool TransA, const Epetra_Vector& x, Epetra_Vector& y, double Alpha) const {
00448   std::cout << "This function Epetra_OskiMatrix::Solve probably works in serial but has not been tested.\n  It will not work in parallel.\n  If you wish to use it feel free to comment out this line and the next return statement.\n  However, correctness and performance are not guaranteed.\n";
00449   return(-1);
00450   Epetra_OskiVector* xCast = NULL;
00451   Epetra_OskiVector* yCast = NULL;
00452   Epetra_OskiVector* tCast = NULL;
00453   bool xCreate = false;
00454   bool yCreate = false;
00455   int ReturnVal;
00456   xCast = dynamic_cast<Epetra_OskiVector*>(const_cast<Epetra_Vector*>(&x));
00457   yCast = dynamic_cast<Epetra_OskiVector*>(&y);
00458   if (xCast == NULL) {
00459     xCast = new Epetra_OskiVector(x);
00460     xCreate = true;
00461   }
00462   if (yCast == NULL) {
00463     yCast = new Epetra_OskiVector(y);
00464     yCreate = true;
00465   }
00466   tCast = new Epetra_OskiVector(x);
00467   if(TransA)
00468     ReturnVal = oski_MatTrisolve(A_tunable_, OP_TRANS, Alpha, (*tCast).Oski_View());
00469   else
00470     ReturnVal = oski_MatTrisolve(A_tunable_, OP_NORMAL, Alpha, (*tCast).Oski_View());
00471   if(ReturnVal)
00472     std::cerr << "OskiVector Solve error\n";
00473   if(xCreate)
00474     delete xCast;
00475   yCast = tCast;
00476   if(yCreate)
00477     delete yCast;
00478   delete tCast;
00479   return ReturnVal;
00480 }
00481 
00482 int Epetra_OskiMatrix::Solve(bool Upper, bool TransA, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00483   int ReturnVal;
00484   ReturnVal = this->Solve(TransA, X, Y, 1.0);
00485   return ReturnVal;
00486 }
00487 
00488 int Epetra_OskiMatrix::Solve(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y, double Alpha) const {
00489   std::cout << "This function Epetra_OskiMatrix::Solve probably works in serial but has not been tested.\n  It will not work in parallel.\n  If you wish to use it feel free to comment out this line and the next return statement.\n  However, correctness and performance are not guaranteed.\n";
00490   return(-1);
00491   Epetra_OskiMultiVector* XCast = NULL;
00492   Epetra_OskiMultiVector* YCast = NULL;
00493   Epetra_OskiMultiVector* TCast = NULL;
00494   bool XCreate = false;
00495   bool YCreate = false;
00496   int ReturnVal;
00497   XCast = dynamic_cast<Epetra_OskiMultiVector*>(const_cast<Epetra_MultiVector*>(&X));
00498   YCast = dynamic_cast<Epetra_OskiMultiVector*>(&Y);
00499   if (XCast == NULL) {
00500     XCast = new Epetra_OskiMultiVector(X);
00501     XCreate = true;
00502   }
00503   if (YCast == NULL) {
00504     YCast = new Epetra_OskiMultiVector(Y);
00505     YCreate = true;
00506   }
00507   TCast = new Epetra_OskiMultiVector(X);
00508   if(TransA)
00509     ReturnVal = oski_MatTrisolve(A_tunable_, OP_TRANS, Alpha, (*TCast).Oski_View());
00510   else
00511     ReturnVal = oski_MatTrisolve(A_tunable_, OP_NORMAL, Alpha, (*TCast).Oski_View());
00512   if(ReturnVal)
00513     std::cerr << "OskiMultiVector Solve error\n";
00514   if(XCreate)
00515     delete XCast;
00516   YCast = TCast;
00517   if(YCreate)
00518     delete YCast;
00519   delete TCast;
00520   return ReturnVal;
00521 }
00522 
00523 int Epetra_OskiMatrix::MatTransMatMultiply(bool ATA, 
00524              const Epetra_Vector& x, 
00525              Epetra_Vector& y, 
00526              Epetra_Vector* t, 
00527              double Alpha, 
00528              double Beta) const {
00529   int ReturnVal;
00530   
00531   if(!Filled())
00532     EPETRA_CHK_ERR(-1); // Matrix must be filled.
00533 
00534   double* xp = (double*) x.Values();
00535   double* xp2 = (double*) x.Values();
00536   double* yp = (double*) y.Values();
00537   double* tp = 0;
00538   if(t != NULL)
00539     tp = (double*) t->Values();
00540 
00541   Epetra_Vector* xcopy = 0;
00542   if (&x==&y && Importer()==0 && Exporter()==0) {
00543     xcopy = new Epetra_Vector(x);
00544     xp = (double *) xcopy->Values();
00545   }
00546   UpdateImportVector(1); // Refresh import and output vectors if needed
00547   UpdateExportVector(1);
00548   
00549 
00550   if(ATA) {
00551     if(Importer() != 0) {
00552       EPETRA_CHK_ERR(ImportVector_->Import(x, *Importer(), Insert));
00553       xp = (double*) ImportVector_->Values();
00554       xp2 = new double[ImportVector_->MyLength()];
00555       for(int i = 0; i < ImportVector_->MyLength(); i++) 
00556         xp2[i] = xp[i];
00557       EPETRA_CHK_ERR(ImportVector_->Import(y, *Importer(), Insert));
00558       yp = (double*) ImportVector_->Values();
00559     }
00560 
00561     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00562     if(Exporter() != 0 && (t) != NULL) {
00563       tp = (double*) ExportVector_->Values();
00564     }
00565 
00566     oski_vecview_t oskiX=0;
00567     oski_vecview_t oskiY=0;
00568     oski_vecview_t oskiT=0;
00569     if(Importer() != 0) 
00570       oskiX = oski_CreateVecView(xp2,ImportVector_->MyLength(),1);
00571     else
00572       oskiX = oski_CreateVecView(xp2,x.MyLength(),1);
00573     if(Importer() != 0) 
00574       oskiY = oski_CreateVecView(yp,ImportVector_->MyLength(),1);
00575     else
00576       oskiY = oski_CreateVecView(yp,y.MyLength(),1);
00577 
00578     if (t != NULL) {
00579       if(Exporter() != 0) 
00580           oskiT = oski_CreateVecView(tp,ExportVector_->MyLength(),1);
00581       else
00582         oskiT = oski_CreateVecView(tp,t->MyLength(),1);
00583 
00584     }
00585     else
00586       oskiT = INVALID_VEC;
00587     ReturnVal = oski_MatTransMatMult(A_tunable_, OP_AT_A, Alpha, oskiX, Beta, oskiY, oskiT);
00588 
00589     if(Importer() != 0) {
00590       y.PutScalar(0.0); // Make sure target is zero
00591       EPETRA_CHK_ERR(y.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
00592     }
00593     // Handle case of rangemap being a local replicated map
00594     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
00595     
00596     if(Exporter() != 0 && (t != NULL)) {
00597       t->PutScalar(0.0); // Make sure target is zero
00598       EPETRA_CHK_ERR(t->Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
00599     }
00600     // Handle case of rangemap being a local replicated map
00601     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(t->Reduce());
00602     UpdateFlops(4 * NumGlobalNonzeros64());
00603 
00604   }
00605   else {
00606     if(this->Comm().NumProc() == 1) {
00607       oski_vecview_t oskiX=0;
00608       oski_vecview_t oskiY=0;
00609       oski_vecview_t oskiT=0;
00610       if (t != NULL)
00611         oskiT = oski_CreateVecView(tp,t->MyLength(),1);
00612       oskiX = oski_CreateVecView(xp,x.MyLength(),1);
00613       oskiY = oski_CreateVecView(yp,y.MyLength(),1);
00614       ReturnVal = oski_MatTransMatMult(A_tunable_, OP_A_AT, Alpha, oskiX, Beta, oskiY, oskiT);
00615       UpdateFlops(4 * NumGlobalNonzeros64());
00616     }
00617     else {
00618 
00619       if(t == NULL) {
00620         Epetra_Vector tempResult(this->DomainMap());
00621         ReturnVal = this->Multiply(false, x, tempResult, 1.0, 0.0);
00622         int ReturnVal2 = this->Multiply(true, tempResult, y, Alpha, Beta);
00623         if(ReturnVal < ReturnVal2)
00624           ReturnVal = ReturnVal2;
00625       }
00626       else {
00627         ReturnVal = this->Multiply(false, x, *t, 1.0, 0.0);
00628         int ReturnVal2 = this->Multiply(true,*t, y, Alpha, Beta);
00629         if(ReturnVal < ReturnVal2)
00630           ReturnVal = ReturnVal2;
00631       }
00632     } 
00633   }
00634   if (xcopy!=0) {
00635     delete xcopy;
00636     EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of x
00637     return(1);
00638   }
00639   return ReturnVal;
00640 }
00641 
00642 int Epetra_OskiMatrix::MatTransMatMultiply(bool ATA, 
00643              const Epetra_MultiVector& X, 
00644              Epetra_MultiVector& Y, 
00645              Epetra_MultiVector* T, 
00646              double Alpha, 
00647                    double Beta) const {
00648   int ReturnVal;
00649   
00650   if(!Filled())
00651     EPETRA_CHK_ERR(-1); // Matrix must be filled.
00652 
00653   int NumVectors = X.NumVectors();
00654   if (NumVectors!=Y.NumVectors()) {
00655     EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
00656   }
00657 
00658   
00659 
00660   double** Xp = (double**) X.Pointers();
00661   double** Xp2 = (double**) X.Pointers();
00662   double** Yp = (double**) Y.Pointers();
00663   double** Tp = 0;
00664   if(T != NULL)
00665     Tp = (double**) T->Pointers();
00666      
00667   int LDX = X.ConstantStride() ? X.Stride() : 0;
00668   int LDY = Y.ConstantStride() ? Y.Stride() : 0;
00669   int LDT = 0;
00670   if(T != NULL)
00671     LDT = T->ConstantStride() ? T->Stride() : 0;
00672 
00673   Epetra_MultiVector* Xcopy = 0;
00674   Epetra_MultiVector* X2 = 0;
00675   if (&X==&Y && Importer()==0 && Exporter()==0) {
00676     Xcopy = new Epetra_MultiVector(X);
00677     Xp = (double **) Xcopy->Pointers();
00678     LDX = Xcopy->ConstantStride() ? Xcopy->Stride() : 0;
00679   }
00680   UpdateImportVector(NumVectors); // Refresh import and output vectors if needed
00681   UpdateExportVector(NumVectors);
00682   
00683 
00684   if(ATA) {
00685     if(Importer() != 0) {
00686       EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
00687       Xp = (double**) ImportVector_->Pointers();
00688       LDX = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
00689 //      need to think about this
00690       X2 = new Epetra_MultiVector(X);
00691       double** Xp2 = (double**) X2->Pointers();
00692       Xp2 = (double **) X2->Pointers();
00693       EPETRA_CHK_ERR(ImportVector_->Import(Y, *Importer(), Insert));
00694       Yp = (double**) ImportVector_->Pointers();
00695       LDY = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
00696     }
00697 
00698     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00699     if(Exporter() != 0 && T != NULL) {
00700       Tp = (double**) ExportVector_->Pointers();
00701       LDT = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
00702     }
00703     oski_vecview_t oskiX=0;
00704     oski_vecview_t oskiY=0;
00705     oski_vecview_t oskiT=0;
00706     if(Importer() != 0) 
00707       oskiX = oski_CreateMultiVecView(*Xp2,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
00708     else
00709       oskiX = oski_CreateMultiVecView(*Xp2,X.MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
00710     if(Importer() != 0) 
00711       oskiY = oski_CreateMultiVecView(*Yp,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
00712     else
00713       oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
00714 
00715     if (T != NULL) {
00716       if(Exporter() != 0) 
00717         oskiT = oski_CreateMultiVecView(*Tp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDT);
00718       else
00719         oskiT = oski_CreateMultiVecView(*Tp,T->MyLength(),NumVectors,LAYOUT_COLMAJ,LDT);
00720 
00721     }
00722     else
00723       oskiT = INVALID_VEC;
00724     ReturnVal = oski_MatTransMatMult(A_tunable_, OP_AT_A, Alpha, oskiX, Beta, oskiY, oskiT);
00725 
00726     if(Importer() != 0) {
00727       Y.PutScalar(0.0); // Make sure target is zero
00728       EPETRA_CHK_ERR(Y.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
00729     }
00730     // Handle case of rangemap being a local replicated map
00731     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
00732     
00733     if(Exporter() != 0 && (T != NULL)) {
00734       T->PutScalar(0.0); // Make sure target is zero
00735       EPETRA_CHK_ERR(T->Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector
00736     }
00737     // Handle case of rangemap being a local replicated map
00738     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(T->Reduce());
00739     UpdateFlops(4 * NumGlobalNonzeros64());
00740 
00741   }
00742   else {
00743     if(this->Comm().NumProc() == 1) {
00744       oski_vecview_t oskiX=0;
00745       oski_vecview_t oskiY=0;
00746       oski_vecview_t oskiT=0;
00747       if (T != NULL)
00748         oskiT = oski_CreateMultiVecView(*Tp,T->MyLength(),NumVectors, LAYOUT_COLMAJ,LDT);
00749       oskiX = oski_CreateMultiVecView(*Xp,X.MyLength(),NumVectors, LAYOUT_COLMAJ,LDX);
00750       oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors, LAYOUT_COLMAJ,LDY);
00751       ReturnVal = oski_MatTransMatMult(A_tunable_, OP_A_AT, Alpha, oskiX, Beta, oskiY, oskiT);
00752       UpdateFlops(4 * NumGlobalNonzeros64() *NumVectors);
00753     }
00754     else {
00755       if(T == NULL) {
00756         Epetra_MultiVector TempResult(this->DomainMap(), NumVectors);
00757         ReturnVal = this->Multiply(false, X, TempResult, 1.0, 0.0);
00758         int ReturnVal2 = this->Multiply(true, TempResult, Y, Alpha, Beta);
00759         if(ReturnVal < ReturnVal2)
00760           ReturnVal = ReturnVal2;
00761       }
00762       else {
00763         ReturnVal = this->Multiply(false, X, *T, 1.0, 0.0);
00764         int ReturnVal2 = this->Multiply(true, *T, Y, Alpha, Beta);
00765         if(ReturnVal < ReturnVal2)
00766           ReturnVal = ReturnVal2;
00767       }
00768     } 
00769   }
00770   if (Xcopy!=0) {
00771     delete Xcopy;
00772     EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of x
00773     return(1);
00774   }
00775   return ReturnVal;
00776 }
00777 
00778 int Epetra_OskiMatrix::MultiplyAndMatTransMultiply(bool TransA, 
00779                const Epetra_Vector& x, 
00780                Epetra_Vector& y, 
00781                const Epetra_Vector& w, 
00782                Epetra_Vector& z, 
00783                double Alpha, 
00784                double Beta, 
00785                double Omega, 
00786                double Zeta) const {
00787 
00788   int ReturnVal;
00789 
00790   if(!Filled())
00791     EPETRA_CHK_ERR(-1); // Matrix must be filled.
00792 
00793   double* xp = (double*) x.Values();
00794   double* xp2 = (double*) x.Values();
00795   double* wp = (double*) w.Values();
00796   double* yp = (double*) y.Values();
00797 //  double* yp2 = (double*) y.Values();
00798   double* zp = (double*) z.Values();
00799   Epetra_MultiVector* yp2 = 0;
00800   Epetra_Vector* xcopy = 0;
00801   if (&x==&y && Importer()==0 && Exporter()==0) {
00802     xcopy = new Epetra_Vector(x);
00803     xp = (double *) xcopy->Values();
00804   }
00805   Epetra_Vector* wcopy = 0;
00806   if (&w==&z && Importer()==0 && Exporter()==0) {
00807     wcopy = new Epetra_Vector(w);
00808     wp = (double *) wcopy->Values();
00809   }
00810   UpdateImportVector(1); // Refresh import and output vectors if needed
00811   UpdateExportVector(1);
00812 
00813   if(TransA) {
00814 
00815     if(Importer() != 0) {
00816       EPETRA_CHK_ERR(ImportVector_->Import(x, *Importer(), Insert));
00817       xp = (double*) ImportVector_->Values();
00818       xp2 = new double[ImportVector_->MyLength()];
00819       for(int i = 0; i < ImportVector_->MyLength(); i++)
00820         xp2[i] = xp[i];
00821       EPETRA_CHK_ERR(ImportVector_->Import(z, *Importer(), Insert));
00822       zp = (double*) ImportVector_->Values();
00823     }
00824     if(Exporter() != 0) {
00825       yp2 = new Epetra_MultiVector(*ExportVector_);
00826       yp = (double*) yp2->Values();
00827       
00828       //for(int i = 0; i < ExportVector_->MyLength(); i++)
00829         //yp2[i] = yp[i];
00830       wp = (double*) ExportVector_->Values();
00831     }
00832 
00833     oski_vecview_t oskiX=0;
00834     oski_vecview_t oskiY=0;
00835     oski_vecview_t oskiW=0;
00836     oski_vecview_t oskiZ=0;
00837     if(Importer() != 0) {
00838       oskiX = oski_CreateVecView(xp2,ImportVector_->MyLength(),1);
00839       oskiZ = oski_CreateVecView(zp,ImportVector_->MyLength(),1);
00840     }
00841     else {
00842       oskiX = oski_CreateVecView(xp2,x.MyLength(),1);
00843       oskiZ = oski_CreateVecView(zp,z.MyLength(),1);
00844     }
00845     if(Exporter() != 0) {
00846       oskiY = oski_CreateVecView(yp,ExportVector_->MyLength(),1);
00847       oskiW = oski_CreateVecView(wp,ExportVector_->MyLength(),1);
00848     }
00849     else {
00850       oskiY = oski_CreateVecView(yp,y.MyLength(),1);
00851       oskiW = oski_CreateVecView(wp,w.MyLength(),1);
00852     }
00853    
00854     ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, oskiX, Beta, oskiY, OP_TRANS, Omega, oskiW, Zeta, oskiZ);
00855     if(Exporter() != 0) {
00856       y.PutScalar(0.0); // Make sure target is zero
00857       EPETRA_CHK_ERR(y.Export(*yp2, *Exporter(), Add)); // Fill y with Values from export vector
00858     }
00859     if(Importer() != 0) {
00860       z.PutScalar(0.0); // Make sure target is zero
00861       EPETRA_CHK_ERR(z.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
00862     }
00863     // Handle case of rangemap being a local replicated map
00864     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
00865     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(z.Reduce());
00866 
00867     UpdateFlops(4 * NumGlobalNonzeros64());
00868   }
00869   //  ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, (*xCast).Oski_View(), Beta, (*yCast).Oski_View(), OP_TRANS, Omega, (*wCast).Oski_View(), Zeta, (*zCast).Oski_View());
00870   else {
00871 
00872    if(Importer() != 0) {
00873       EPETRA_CHK_ERR(ImportVector_->Import(x, *Importer(), Insert));
00874       xp = (double*) ImportVector_->Values();
00875       xp2 = new double[ImportVector_->MyLength()];
00876       for(int i = 0; i < ImportVector_->MyLength(); i++)
00877         xp2[i] = xp[i];
00878       EPETRA_CHK_ERR(ImportVector_->Import(w, *Importer(), Insert));
00879       wp = (double*) ImportVector_->Values();
00880     }
00881     if(Exporter() != 0) {
00882       yp2 = new Epetra_MultiVector(*ExportVector_);
00883       yp = (double*) yp2->Values();
00884       
00885       //for(int i = 0; i < ExportVector_->MyLength(); i++)
00886         //yp2[i] = yp[i];
00887       zp = (double*) ExportVector_->Values();
00888     }
00889 
00890     oski_vecview_t oskiX=0;
00891     oski_vecview_t oskiY=0;
00892     oski_vecview_t oskiW=0;
00893     oski_vecview_t oskiZ=0;
00894     if(Importer() != 0) {
00895       oskiX = oski_CreateVecView(xp2,ImportVector_->MyLength(),1);
00896       oskiW = oski_CreateVecView(wp,ImportVector_->MyLength(),1);
00897     }
00898     else {
00899       oskiX = oski_CreateVecView(xp2,x.MyLength(),1);
00900       oskiW = oski_CreateVecView(wp,w.MyLength(),1);
00901     }
00902     if(Exporter() != 0) {
00903       oskiY = oski_CreateVecView(yp,ExportVector_->MyLength(),1);
00904       oskiZ = oski_CreateVecView(zp,ExportVector_->MyLength(),1);
00905     }
00906     else {
00907       oskiY = oski_CreateVecView(yp,y.MyLength(),1);
00908       oskiZ = oski_CreateVecView(zp,z.MyLength(),1);
00909     }
00910    
00911     ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, oskiX, Beta, oskiY, OP_NORMAL, Omega, oskiW, Zeta, oskiZ);
00912     if(Exporter() != 0) {
00913       y.PutScalar(0.0); // Make sure target is zero
00914       EPETRA_CHK_ERR(y.Export(*yp2, *Exporter(), Add)); // Fill y with Values from export vector
00915       z.PutScalar(0.0); // Make sure target is zero
00916       EPETRA_CHK_ERR(z.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector*/
00917     }
00918     // Handle case of rangemap being a local replicated map
00919     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(y.Reduce());
00920     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(z.Reduce());
00921 
00922     UpdateFlops(4 * NumGlobalNonzeros64());
00923 
00924   }
00925   if(ReturnVal)
00926     std::cerr << "OskiVector multiply error\n";
00927   return ReturnVal;
00928 }
00929 
00930 int Epetra_OskiMatrix::MultiplyAndMatTransMultiply(bool TransA, 
00931                const Epetra_MultiVector& X, 
00932                Epetra_MultiVector& Y, 
00933                const Epetra_MultiVector& W, 
00934                Epetra_MultiVector& Z, 
00935                double Alpha, 
00936                double Beta, 
00937                double Omega, 
00938                double Zeta) const {
00939   int ReturnVal;
00940   if(!Filled())
00941     EPETRA_CHK_ERR(-1); // Matrix must be filled.
00942   int NumVectors = X.NumVectors();
00943   if (NumVectors!=Y.NumVectors()) {
00944     EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
00945   }
00946 
00947 
00948   double** Xp = (double**) X.Pointers();
00949   double** Xp2 = (double**) X.Pointers();
00950   double** Wp = (double**) W.Pointers();
00951   double** Yp = (double**) Y.Pointers();
00952   double** Zp = (double**) Z.Pointers();
00953   int LDX = X.ConstantStride() ? X.Stride() : 0;
00954   int LDY = Y.ConstantStride() ? Y.Stride() : 0;
00955   int LDW = W.ConstantStride() ? W.Stride() : 0;
00956   int LDZ = Z.ConstantStride() ? Z.Stride() : 0;
00957 
00958   Epetra_MultiVector* Yp2 = 0;
00959   Epetra_MultiVector* X2 = 0;
00960   Epetra_MultiVector* Xcopy = 0;
00961   if (&X==&Y && Importer()==0 && Exporter()==0) {
00962     Xcopy = new Epetra_MultiVector(X);
00963     Xp = (double **) Xcopy->Pointers();
00964     LDX = Xcopy->ConstantStride() ? Xcopy->Stride() : 0;
00965   }
00966   Epetra_MultiVector* Wcopy = 0;
00967   if (&W==&Z && Importer()==0 && Exporter()==0) {
00968     Wcopy = new Epetra_MultiVector(W);
00969     Wp = (double **) Wcopy->Values();
00970     LDW = Wcopy->ConstantStride() ? Wcopy->Stride() : 0;
00971   }
00972   UpdateImportVector(NumVectors); // Refresh import and output vectors if needed
00973   UpdateExportVector(NumVectors);
00974 
00975   if(TransA) {
00976 
00977     if(Importer() != 0) {
00978       EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
00979       Xp = (double**) ImportVector_->Pointers();
00980       LDX = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
00981       X2 = new Epetra_MultiVector(X);
00982       double** Xp2 = (double**) X2->Pointers();
00983       Xp2 = (double **) X2->Pointers();
00984       EPETRA_CHK_ERR(ImportVector_->Import(Y, *Importer(), Insert));
00985       Zp = (double**) ImportVector_->Pointers();
00986       LDZ = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
00987     }
00988     if(Exporter() != 0) {
00989       Yp2 = new Epetra_MultiVector(*ExportVector_);
00990       Yp = (double**) Yp2->Pointers();
00991       Wp = (double**) ExportVector_->Pointers();
00992     }
00993 
00994     oski_vecview_t oskiX=0;
00995     oski_vecview_t oskiY=0;
00996     oski_vecview_t oskiW=0;
00997     oski_vecview_t oskiZ=0;
00998     if(Importer() != 0) {
00999       oskiX = oski_CreateMultiVecView(*Xp2,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
01000       oskiZ = oski_CreateMultiVecView(*Zp,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDZ);
01001     }
01002     else {
01003       oskiX = oski_CreateMultiVecView(*Xp2,X.MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
01004       oskiZ = oski_CreateMultiVecView(*Zp,Z.MyLength(),NumVectors,LAYOUT_COLMAJ,LDZ);
01005     }
01006     if(Exporter() != 0) {
01007       oskiY = oski_CreateMultiVecView(*Yp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
01008       oskiW = oski_CreateMultiVecView(*Wp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDW);
01009     }
01010     else {
01011       oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
01012       oskiW = oski_CreateMultiVecView(*Wp,W.MyLength(),NumVectors,LAYOUT_COLMAJ,LDW);
01013     }
01014    
01015     ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, oskiX, Beta, oskiY, OP_TRANS, Omega, oskiW, Zeta, oskiZ);
01016     if(Exporter() != 0) {
01017       Y.PutScalar(0.0); // Make sure target is zero
01018       EPETRA_CHK_ERR(Y.Export(*Yp2, *Exporter(), Add)); // Fill y with Values from export vector
01019     }
01020     if(Importer() != 0) {
01021       Z.PutScalar(0.0); // Make sure target is zero
01022       EPETRA_CHK_ERR(Z.Export(*ImportVector_, *Importer(), Add)); // Fill y with Values from export vector
01023     }
01024     // Handle case of rangemap being a local replicated map
01025     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
01026     if (!Graph().DomainMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Z.Reduce());
01027 
01028     UpdateFlops(4 * NumGlobalNonzeros64());
01029   }
01030   //  ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, (*xCast).Oski_View(), Beta, (*yCast).Oski_View(), OP_TRANS, Omega, (*wCast).Oski_View(), Zeta, (*zCast).Oski_View());
01031   else {
01032 
01033    if(Importer() != 0) {
01034       EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
01035       Xp = (double**) ImportVector_->Pointers();
01036       LDX = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
01037       X2 = new Epetra_MultiVector(X);
01038       Xp2 = (double**) X2->Pointers();
01039       EPETRA_CHK_ERR(ImportVector_->Import(W, *Importer(), Insert));
01040       Wp = (double**) ImportVector_->Pointers();
01041       LDW = ImportVector_->ConstantStride() ? ImportVector_->Stride() : 0;
01042     }
01043     if(Exporter() != 0) {
01044       Yp2 = new Epetra_MultiVector(*ExportVector_);
01045       Yp = (double**) Yp2->Pointers();
01046       Zp = (double**) ExportVector_->Pointers();
01047     }
01048 
01049     oski_vecview_t oskiX=0;
01050     oski_vecview_t oskiY=0;
01051     oski_vecview_t oskiW=0;
01052     oski_vecview_t oskiZ=0;
01053     if(Importer() != 0) {
01054       oskiX = oski_CreateMultiVecView(*Xp2,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
01055       oskiW = oski_CreateMultiVecView(*Wp,ImportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDW);
01056     }
01057     else {
01058       oskiX = oski_CreateMultiVecView(*Xp2,X.MyLength(),NumVectors,LAYOUT_COLMAJ,LDX);
01059       oskiW = oski_CreateMultiVecView(*Wp,W.MyLength(),NumVectors,LAYOUT_COLMAJ,LDW);
01060     }
01061     if(Exporter() != 0) {
01062       oskiY = oski_CreateMultiVecView(*Yp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
01063       oskiZ = oski_CreateMultiVecView(*Zp,ExportVector_->MyLength(),NumVectors,LAYOUT_COLMAJ,LDZ);
01064     }
01065     else {
01066       oskiY = oski_CreateMultiVecView(*Yp,Y.MyLength(),NumVectors,LAYOUT_COLMAJ,LDY);
01067       oskiZ = oski_CreateMultiVecView(*Zp,Z.MyLength(),NumVectors,LAYOUT_COLMAJ,LDZ);
01068     }
01069    
01070     ReturnVal = oski_MatMultAndMatTransMult(A_tunable_, Alpha, oskiX, Beta, oskiY, OP_NORMAL, Omega, oskiW, Zeta, oskiZ);
01071     if(Exporter() != 0) {
01072       Y.PutScalar(0.0); // Make sure target is zero
01073       EPETRA_CHK_ERR(Y.Export(*Yp2, *Exporter(), Add)); // Fill y with Values from export vector
01074       Z.PutScalar(0.0); // Make sure target is zero
01075       EPETRA_CHK_ERR(Z.Export(*ExportVector_, *Exporter(), Add)); // Fill y with Values from export vector*/
01076     }
01077     // Handle case of rangemap being a local replicated map
01078     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
01079     if (!Graph().RangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Z.Reduce());
01080 
01081     UpdateFlops(4 * NumGlobalNonzeros64());
01082 
01083   }
01084   return ReturnVal;
01085 }
01086 
01087 int Epetra_OskiMatrix::MatPowMultiply(bool TransA, 
01088               const Epetra_Vector& x, 
01089               Epetra_Vector& y, 
01090               Epetra_MultiVector& T, 
01091               int Power,
01092               double Alpha, 
01093               double Beta) const {
01094   //The code has not been tested.  It should work in serial but not in parallel.
01095   std::cerr << "MatPowMultiply not implemented in oski-1.01h release.\n";
01096   return -1;
01097 
01098   int ReturnVal; 
01099 
01100   if(!Filled())
01101     EPETRA_CHK_ERR(-1); // Matrix must be filled.
01102 
01103   double* xp = (double*) x.Values();
01104   double* yp = (double*) y.Values();
01105   double** Tp = (double**) T.Pointers();
01106 
01107   Epetra_MultiVector *Tptr;
01108 
01109   int LDT = T.ConstantStride() ? T.Stride() : 0;
01110 
01111 
01112   if(this->Comm().NumProc() == 1) {
01113     oski_vecview_t oskiX=0;
01114     oski_vecview_t oskiY=0;
01115     oski_vecview_t oskiT=0;
01116     if (&T != NULL)
01117       oskiT = oski_CreateMultiVecView(*Tp,T.MyLength(),Power,LAYOUT_COLMAJ,LDT);
01118     oskiX = oski_CreateVecView(xp,x.MyLength(),1);
01119     oskiY = oski_CreateVecView(yp,y.MyLength(),1);
01120 
01121     if(TransA)
01122       ReturnVal = oski_MatPowMult(A_tunable_, OP_TRANS, Power, Alpha, oskiX, Beta, oskiY, oskiT);
01123     else
01124       ReturnVal = oski_MatPowMult(A_tunable_, OP_NORMAL, Power, Alpha, oskiX, Beta, oskiY, oskiT);
01125     if(ReturnVal)
01126       std::cerr << "OskiVector matpow multiply error\n";
01127   }
01128   else {
01129     if(&T == NULL) 
01130       Tptr = new Epetra_MultiVector(x.Map(), Power);
01131     else
01132       Tptr = &T;
01133     if(TransA) {
01134       ReturnVal = this->Multiply(true, x, Tptr[1], 1.0, 0.0);
01135       for(int i = 1; i < Power-1; i++)
01136         ReturnVal = this->Multiply(true, Tptr[i], Tptr[i+1], 1.0, 0.0);
01137       ReturnVal = this->Multiply(true, Tptr[Power-2], y, Alpha, Beta);
01138     }
01139     else {
01140       ReturnVal = this->Multiply(false, x, Tptr[1], 1.0, 0.0);
01141       for(int i = 1; i < Power-1; i++)
01142         ReturnVal = this->Multiply(false, Tptr[i], Tptr[i+1], 1.0, 0.0);
01143       ReturnVal = this->Multiply(false, Tptr[Power-2], y, Alpha, Beta);
01144     }
01145     if(ReturnVal)
01146       std::cerr << "OskiVector matpow multiply error\n";
01147     if(&T == NULL)
01148       delete(Tptr);
01149   }    
01150   UpdateFlops(2 * Power * NumGlobalNonzeros64());
01151   return ReturnVal;
01152 }
01153 
01154 int Epetra_OskiMatrix::MatPowMultiply(bool TransA, 
01155               const Epetra_Vector& x, 
01156               Epetra_Vector& y, 
01157               int Power,
01158               double Alpha, 
01159               double Beta) const {
01160 
01161   //The code has not been tested.  It should work in serial but not in parallel.
01162   std::cerr << "MatPowMultiply not implemented in oski-1.01h release.\n";
01163   return -1;
01164   std::cerr << "mult\n";
01165   Epetra_OskiVector* xCast = NULL;
01166   Epetra_OskiVector* yCast = NULL;
01167   bool xCreate = false;
01168   bool yCreate = false;
01169   int ReturnVal;
01170   xCast = dynamic_cast<Epetra_OskiVector*>(const_cast <Epetra_Vector*>(&x));
01171   yCast = dynamic_cast<Epetra_OskiVector*>(&y);
01172   if (xCast == NULL) {
01173     xCast = new Epetra_OskiVector(x);
01174     xCreate = true;
01175   }
01176   if (yCast == NULL) {
01177     yCast = new Epetra_OskiVector(y);
01178     yCreate = true;
01179   }
01180   std::cerr << "mult\n";
01181   if(TransA)
01182     ReturnVal = oski_MatPowMult(A_tunable_, OP_TRANS, Power, Alpha, (*xCast).Oski_View(), Beta, (*yCast).Oski_View(), NULL);
01183   else
01184     ReturnVal = oski_MatPowMult(A_tunable_, OP_NORMAL, Power, Alpha, (*xCast).Oski_View(), Beta, (*yCast).Oski_View(), NULL);
01185   std::cerr << "done\n";
01186   if(ReturnVal)
01187     std::cerr << "OskiVector matpow multiply error\n";
01188   if(xCreate)
01189     delete xCast;
01190   if(yCreate)
01191     delete yCast;
01192   return ReturnVal;
01193 }
01194 
01195 int Epetra_OskiMatrix::SetHint(const Teuchos::ParameterList& List) {
01196   int* ArgArray = NULL;
01197   int Diags;
01198   int Blocks;
01199   char Number[10];
01200   char Row[20];
01201   char Col[20];
01202   char Diag[20];
01203   int ReturnVal = 0;
01204   if(List.isParameter("randompattern"))
01205     if(Teuchos::getParameter<bool>(List, "randompattern"))
01206       if(ReturnVal = oski_SetHint(A_tunable_, HINT_RANDOM_PATTERN))
01207         std::cerr << "Error setting hint random pattern.\n";
01208   if(List.isParameter("correlatedpattern"))
01209     if(Teuchos::getParameter<bool>(List, "correlatedpattern"))
01210       if(ReturnVal = oski_SetHint(A_tunable_, HINT_CORREL_PATTERN))
01211         std::cerr << "Error setting hint correlated pattern.\n";
01212   if(List.isParameter("symmetricpattern"))
01213     if(Teuchos::getParameter<bool>(List, "symmetricpattern"))
01214       if(ReturnVal = oski_SetHint(A_tunable_, HINT_SYMM_PATTERN))
01215         std::cerr << "Error setting hint symmetric pattern.\n";
01216   if(List.isParameter("nonsymmetricpattern"))
01217     if(Teuchos::getParameter<bool>(List, "nonsymmetricpattern"))
01218       if(ReturnVal = oski_SetHint(A_tunable_, HINT_NONSYMM_PATTERN))
01219         std::cerr << "Error setting hint nonsymmetric pattern.\n";
01220   if(List.isParameter("alignedblocks"))
01221     if(Teuchos::getParameter<bool>(List, "alignedblocks"))
01222     {
01223       if(ReturnVal = oski_SetHint(A_tunable_, HINT_ALIGNED_BLOCKS))
01224         std::cerr << "Error setting hint aligned blocks.\n";
01225     }
01226   if(List.isParameter("unalignedblocks"))
01227     if(Teuchos::getParameter<bool>(List, "unalignedblocks"))
01228       if(ReturnVal = oski_SetHint(A_tunable_, HINT_UNALIGNED_BLOCKS))
01229         std::cerr << "Error setting hint unaligned blocks.\n";
01230   if(List.isParameter("nodiags"))
01231     if(Teuchos::getParameter<bool>(List, "nodiags"))
01232       if(ReturnVal = oski_SetHint(A_tunable_, HINT_NO_DIAGS))
01233         std::cerr << "Error setting hint no diags.\n";
01234   if(List.isParameter("noblocks"))
01235     if(Teuchos::getParameter<bool>(List, "noblocks"))
01236       if(ReturnVal = oski_SetHint(A_tunable_, HINT_NO_BLOCKS))
01237         std::cerr << "Error setting hint no blocks.\n";
01238   if(List.isParameter("singleblocksize"))
01239     if(Teuchos::getParameter<bool>(List, "singleblocksize")) {
01240       ArgArray = new int[2];
01241       if(List.isParameter("row"))
01242       {
01243         ArgArray[0] = Teuchos::getParameter<int>(List, "row");
01244         if(List.isParameter("col"))
01245           ArgArray[1] = Teuchos::getParameter<int>(List, "col");
01246         if(ReturnVal = oski_SetHint(A_tunable_, HINT_SINGLE_BLOCKSIZE, ArgArray[0], ArgArray[1]))
01247           std::cerr << "Error setting hint single block size.\n";
01248       }
01249       else
01250         if(ReturnVal = oski_SetHint(A_tunable_, HINT_SINGLE_BLOCKSIZE))
01251           std::cerr << "Error setting hint single block size.\n";
01252       delete [] ArgArray;
01253       ArgArray = NULL;
01254     }
01255   if(List.isParameter("multipleblocksize"))
01256     if(Teuchos::getParameter<bool>(List, "multipleblocksize"))
01257       if(List.isParameter("blocks")) {
01258         Blocks = Teuchos::getParameter<int>(List, "blocks");
01259         ArgArray = new int[2*Blocks+1];
01260   ArgArray[0] = Blocks;
01261         for(int i = 0; i < Blocks; i++) {
01262     sprintf(Number, "%d", i+1);
01263     strcpy(Row, "row");
01264     strcpy(Col, "col");
01265     strcat(Row, Number);
01266     strcat(Col, Number);
01267           if(List.isParameter(Row))
01268             ArgArray[i*2 + 1] = Teuchos::getParameter<int>(List, Row);
01269           if(List.isParameter(Col))
01270             ArgArray[i*2 + 2] = Teuchos::getParameter<int>(List, Col);
01271         }
01272         switch(Blocks) {
01273           case 1 :  if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES, ArgArray[0], ArgArray[1], ArgArray[2]))
01274                   std::cerr << "Error setting hint multiple blocks.\n";
01275         break;
01276           case 2 :  if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4]))
01277                   std::cerr << "Error setting hint multiple blocks.\n";
01278         break;
01279           case 3 :  if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4], ArgArray[5], ArgArray[6]))
01280                   std::cerr << "Error setting hint multiple blocks.\n";
01281         break;
01282           case 4 :  if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4], ArgArray[5], ArgArray[6], ArgArray[7], ArgArray[8]))
01283                   std::cerr << "Error setting hint multiple blocks.\n";
01284         break;
01285           case 5 :  if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4], ArgArray[5], ArgArray[6], ArgArray[7], ArgArray[8], ArgArray[9], ArgArray[10]))
01286                   std::cerr << "Error setting hint multiple blocks.\n";
01287         break;
01288           default : if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES))
01289                       std::cerr << "Error setting hint multiple blocks.\n";
01290                     break;
01291         }
01292         delete [] ArgArray;
01293         ArgArray = NULL;
01294       }
01295       else
01296         if(ReturnVal = oski_SetHint(A_tunable_, HINT_MULTIPLE_BLOCKSIZES))
01297           std::cerr << "Error setting hint multiple blocks.\n";
01298   if(List.isParameter("diags"))
01299     if(Teuchos::getParameter<bool>(List, "diags"))
01300       if(List.isParameter("numdiags")) {
01301         Diags = Teuchos::getParameter<int>(List, "numdiags");
01302         ArgArray = new int[Diags+1];
01303   ArgArray[0] = Diags;
01304         for(int i = 0; i < Diags; i++) {
01305     sprintf(Number, "%d", i + 1);
01306     strcpy(Diag, "diag");
01307     strcat(Diag, Number);
01308           if(List.isParameter(Diag))
01309             ArgArray[i + 1] = Teuchos::getParameter<int>(List, Diag);
01310         }
01311         switch(Diags) {
01312           case 1 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0], ArgArray[1]))
01313                      std::cerr << "Error setting hint diags\n";
01314                    break;
01315           case 2 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0], ArgArray[1], ArgArray[2]))
01316                      std::cerr << "Error setting hint diags\n";
01317                    break;
01318           case 3 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3]))
01319                      std::cerr << "Error setting hint diags\n";
01320                    break;
01321           case 4 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4]))
01322                      std::cerr << "Error setting hint diags\n";
01323                    break;
01324           case 5 : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0], ArgArray[1], ArgArray[2], ArgArray[3], ArgArray[4], ArgArray[5]))
01325                      std::cerr << "Error setting hint diags\n";
01326                    break;
01327           default : if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS, ArgArray[0]))
01328                       std::cerr << "Error setting hint diags\n";
01329                     break;
01330         }
01331         delete [] ArgArray;
01332       }
01333       else
01334       {
01335         if(ReturnVal = oski_SetHint(A_tunable_, HINT_DIAGS))
01336           std::cerr << "Error setting hint digs.\n";
01337       }
01338   return ReturnVal;
01339 }
01340 
01341 int Epetra_OskiMatrix::SetHintMultiply(bool TransA, 
01342                double Alpha, 
01343                const Epetra_OskiMultiVector& InVec, 
01344                double Beta, 
01345                const Epetra_OskiMultiVector& OutVec, 
01346                int NumCalls, 
01347                const Teuchos::ParameterList& List) {
01348   int ReturnVal;
01349   oski_vecview_t InView = NULL;
01350   oski_vecview_t OutView = NULL;
01351   InView = InVec.Oski_View();
01352   OutView = OutVec.Oski_View();
01353   if(List.isParameter("tune"))
01354     if(Teuchos::getParameter<bool>(List, "tune"))
01355       NumCalls = ALWAYS_TUNE;
01356   if(List.isParameter("tuneaggressive"))
01357     if(Teuchos::getParameter<bool>(List, "tuneaggressive"))
01358       NumCalls = ALWAYS_TUNE_AGGRESSIVELY;
01359   if(List.isParameter("symminvec"))
01360     if(Teuchos::getParameter<bool>(List, "symminvec"))
01361       InView = SYMBOLIC_VEC;
01362   if(List.isParameter("symminmultivec"))
01363     if(Teuchos::getParameter<bool>(List, "symminmultivec"))
01364       InView = SYMBOLIC_MULTIVEC;
01365   if(List.isParameter("symmoutvec"))
01366     if(Teuchos::getParameter<bool>(List, "symmoutvec"))
01367       OutView = SYMBOLIC_VEC;
01368   if(List.isParameter("symmoutmultivec"))
01369     if(Teuchos::getParameter<bool>(List, "symmoutmultivec"))
01370       OutView = SYMBOLIC_MULTIVEC;
01371   if(this->Comm().NumProc() > 1) {
01372     if(InVec.NumVectors() == 1)
01373       InView = SYMBOLIC_VEC;
01374     else
01375       InView = SYMBOLIC_MULTIVEC;
01376     if(OutVec.NumVectors() == 1)
01377       OutView = SYMBOLIC_VEC;
01378     else
01379       OutView = SYMBOLIC_MULTIVEC;
01380   }
01381   if(TransA)
01382     ReturnVal = oski_SetHintMatMult(A_tunable_, OP_TRANS, Alpha, InView, Beta, OutView, NumCalls);
01383   else
01384     ReturnVal = oski_SetHintMatMult(A_tunable_, OP_NORMAL, Alpha, InView, Beta, OutView, NumCalls);
01385   if(ReturnVal)
01386     std::cerr << "Set hint multivector multiply error\n";
01387   return ReturnVal;
01388 }
01389 
01390 int Epetra_OskiMatrix::SetHintSolve(bool TransA,
01391             double Alpha, 
01392             const Epetra_OskiMultiVector& Vector, 
01393             int NumCalls, 
01394             const Teuchos::ParameterList& List) {
01395   int ReturnVal;
01396   oski_vecview_t VecView = NULL;
01397   VecView = Vector.Oski_View();
01398   if(List.isParameter("tune"))
01399     if(Teuchos::getParameter<bool>(List, "tune"))
01400       NumCalls = ALWAYS_TUNE;
01401   if(List.isParameter("tuneaggressive"))
01402     if(Teuchos::getParameter<bool>(List, "tuneaggressive"))
01403       NumCalls = ALWAYS_TUNE_AGGRESSIVELY;
01404   if(List.isParameter("symmvec"))
01405     if(Teuchos::getParameter<bool>(List, "symmvec"))
01406       VecView = SYMBOLIC_VEC;
01407   if(List.isParameter("symmmultivec"))
01408     if(Teuchos::getParameter<bool>(List, "symmmultivec"))
01409       VecView = SYMBOLIC_MULTIVEC;
01410   if(this->Comm().NumProc() > 1) {
01411     if(Vector.NumVectors() == 1)
01412       VecView = SYMBOLIC_VEC;
01413     else
01414       VecView = SYMBOLIC_MULTIVEC;
01415   }
01416   if(TransA)
01417     ReturnVal = oski_SetHintMatTrisolve(A_tunable_, OP_TRANS, Alpha, VecView, NumCalls);
01418   else
01419     ReturnVal = oski_SetHintMatTrisolve(A_tunable_, OP_NORMAL, Alpha, VecView, NumCalls);
01420   if(ReturnVal)
01421     std::cerr << "Set hint solve error\n";
01422   return ReturnVal;
01423 }
01424 
01425 int Epetra_OskiMatrix::SetHintMatTransMatMultiply (bool ATA, 
01426                double Alpha, 
01427                const Epetra_OskiMultiVector& InVec, 
01428                double Beta, 
01429                const Epetra_OskiMultiVector& OutVec, 
01430                const Epetra_OskiMultiVector& Intermediate, 
01431                int NumCalls, 
01432                const Teuchos::ParameterList& List) {
01433   int ReturnVal;
01434   oski_vecview_t InView = NULL;
01435   oski_vecview_t OutView = NULL;
01436   oski_vecview_t IntermediateView = NULL;
01437   InView = InVec.Oski_View();
01438   OutView = OutVec.Oski_View();
01439   if(&Intermediate != NULL)
01440     IntermediateView = Intermediate.Oski_View();
01441   if(List.isParameter("tune"))
01442     if(Teuchos::getParameter<bool>(List, "tune"))
01443       NumCalls = ALWAYS_TUNE;
01444   if(List.isParameter("tuneaggressive"))
01445     if(Teuchos::getParameter<bool>(List, "tuneaggressive"))
01446       NumCalls = ALWAYS_TUNE_AGGRESSIVELY;
01447   if(List.isParameter("symminvec"))
01448     if(Teuchos::getParameter<bool>(List, "symminvec"))
01449       InView = SYMBOLIC_VEC;
01450   if(List.isParameter("symminmultivec"))
01451     if(Teuchos::getParameter<bool>(List, "symminmultivec"))
01452       InView = SYMBOLIC_MULTIVEC;
01453   if(List.isParameter("symmoutvec"))
01454     if(Teuchos::getParameter<bool>(List, "symmoutvec"))
01455       OutView = SYMBOLIC_VEC;
01456   if(List.isParameter("symmoutmultivec"))
01457     if(Teuchos::getParameter<bool>(List, "symmoutmultivec"))
01458       OutView = SYMBOLIC_MULTIVEC;
01459   if(List.isParameter("symmintervec"))
01460     if(Teuchos::getParameter<bool>(List, "symmintervec"))
01461       IntermediateView = SYMBOLIC_VEC;
01462   if(List.isParameter("symmintermultivec"))
01463     if(Teuchos::getParameter<bool>(List, "symmintermultivec"))
01464       IntermediateView = SYMBOLIC_MULTIVEC;
01465   if(List.isParameter("invalidinter"))
01466     if(Teuchos::getParameter<bool>(List, "invalidinter"))
01467       IntermediateView = NULL;
01468   if(this->Comm().NumProc() > 1) {
01469     if(InVec.NumVectors() == 1)
01470       InView = SYMBOLIC_VEC;
01471     else
01472       InView = SYMBOLIC_MULTIVEC;
01473     if(OutVec.NumVectors() == 1)
01474       OutView = SYMBOLIC_VEC;
01475     else
01476       OutView = SYMBOLIC_MULTIVEC;
01477     if(Intermediate.NumVectors() == 1)
01478       IntermediateView = SYMBOLIC_VEC;
01479     else
01480       IntermediateView = SYMBOLIC_MULTIVEC;
01481   }
01482   if(ATA)
01483     ReturnVal = oski_SetHintMatTransMatMult(A_tunable_, OP_AT_A, Alpha, InView, Beta, OutView, IntermediateView, NumCalls);
01484   else
01485     ReturnVal = oski_SetHintMatTransMatMult(A_tunable_, OP_A_AT, Alpha, InView, Beta, OutView, IntermediateView, NumCalls);
01486   if(ReturnVal)
01487     std::cerr << "Set hint mattransmat multiply error\n";
01488   return ReturnVal;
01489 }
01490 
01491 int Epetra_OskiMatrix::SetHintMultiplyAndMatTransMultiply(bool TransA, 
01492                 double Alpha, 
01493                 const Epetra_OskiMultiVector& InVec, 
01494                 double Beta, 
01495                 const Epetra_OskiMultiVector& OutVec, 
01496                 double Omega, 
01497                 const Epetra_OskiMultiVector& InVec2, 
01498                 double Zeta, 
01499                 const Epetra_OskiMultiVector& OutVec2, 
01500                 int NumCalls, 
01501                 const Teuchos::ParameterList& List) {
01502   int ReturnVal;
01503   oski_vecview_t InView = NULL;
01504   oski_vecview_t OutView = NULL;
01505   oski_vecview_t InView2 = NULL;
01506   oski_vecview_t OutView2 = NULL;
01507   InView = InVec.Oski_View();
01508   OutView = OutVec.Oski_View();
01509   InView2 = InVec2.Oski_View();
01510   OutView2 = OutVec2.Oski_View();
01511   if(List.isParameter("tune"))
01512     if(Teuchos::getParameter<bool>(List, "tune"))
01513       NumCalls = ALWAYS_TUNE;
01514   if(List.isParameter("tuneaggressive"))
01515     if(Teuchos::getParameter<bool>(List, "tuneaggressive"))
01516       NumCalls = ALWAYS_TUNE_AGGRESSIVELY;
01517   if(List.isParameter("symminvec"))
01518     if(Teuchos::getParameter<bool>(List, "symminvec"))
01519       InView = SYMBOLIC_VEC;
01520   if(List.isParameter("symminmultivec"))
01521     if(Teuchos::getParameter<bool>(List, "symminmultivec"))
01522       InView = SYMBOLIC_MULTIVEC;
01523   if(List.isParameter("symmoutvec"))
01524     if(Teuchos::getParameter<bool>(List, "symmoutvec"))
01525       OutView = SYMBOLIC_VEC;
01526   if(List.isParameter("symmoutmultivec"))
01527     if(Teuchos::getParameter<bool>(List, "symmoutmultivec"))
01528       OutView = SYMBOLIC_MULTIVEC;
01529   if(List.isParameter("symminvec2"))
01530     if(Teuchos::getParameter<bool>(List, "symminvec2"))
01531       InView2 = SYMBOLIC_VEC;
01532   if(List.isParameter("symminmultivec2"))
01533     if(Teuchos::getParameter<bool>(List, "symminmultivec2"))
01534       InView2 = SYMBOLIC_MULTIVEC;
01535   if(List.isParameter("symmoutvec2"))
01536     if(Teuchos::getParameter<bool>(List, "symmoutvec2"))
01537       OutView2 = SYMBOLIC_VEC;
01538   if(List.isParameter("symmoutmultivec2"))
01539     if(Teuchos::getParameter<bool>(List, "symmoutmultivec2"))
01540       OutView2 = SYMBOLIC_MULTIVEC;
01541   if(this->Comm().NumProc() > 1) {
01542     if(InVec.NumVectors() == 1)
01543       InView = SYMBOLIC_VEC;
01544     else
01545       InView = SYMBOLIC_MULTIVEC;
01546     if(OutVec.NumVectors() == 1)
01547       OutView = SYMBOLIC_VEC;
01548     else
01549       OutView = SYMBOLIC_MULTIVEC;
01550     if(InVec2.NumVectors() == 1)
01551       InView2 = SYMBOLIC_VEC;
01552     else
01553       InView2 = SYMBOLIC_MULTIVEC;
01554     if(OutVec2.NumVectors() == 1)
01555       OutView2 = SYMBOLIC_VEC;
01556     else
01557       OutView2 = SYMBOLIC_MULTIVEC;
01558   }
01559   if(TransA)
01560     ReturnVal = oski_SetHintMatMultAndMatTransMult(A_tunable_, Alpha, InView, Beta, OutView, OP_TRANS, Omega, InView2, Zeta, OutView2, NumCalls);
01561   else
01562     ReturnVal = oski_SetHintMatMultAndMatTransMult(A_tunable_, Alpha, InView, Beta, OutView, OP_NORMAL, Omega, InView2, Zeta, OutView2, NumCalls);
01563   if(ReturnVal)
01564     std::cerr << "Set hint multivector multiply and mattransmultiply error\n";
01565   return ReturnVal;
01566 }
01567 
01568 int Epetra_OskiMatrix::SetHintPowMultiply(bool TransA, 
01569             double Alpha, 
01570             const Epetra_OskiMultiVector& InVec, 
01571             double Beta, 
01572             const Epetra_OskiMultiVector& OutVec, 
01573             const Epetra_OskiMultiVector& Intermediate, 
01574             int Power, 
01575             int NumCalls, 
01576             const Teuchos::ParameterList& List) {
01577   int ReturnVal;
01578   oski_vecview_t InView = NULL;
01579   oski_vecview_t OutView = NULL;
01580   oski_vecview_t IntermediateView = NULL;
01581   InView = InVec.Oski_View();
01582   OutView = OutVec.Oski_View();
01583   if(&Intermediate != NULL)
01584     IntermediateView = Intermediate.Oski_View();
01585   if(List.isParameter("tune"))
01586     if(Teuchos::getParameter<bool>(List, "tune"))
01587       NumCalls = ALWAYS_TUNE;
01588   if(List.isParameter("tuneaggressive"))
01589     if(Teuchos::getParameter<bool>(List, "tuneaggressive"))
01590       NumCalls = ALWAYS_TUNE_AGGRESSIVELY;
01591   if(List.isParameter("symminvec"))
01592     if(Teuchos::getParameter<bool>(List, "symminvec"))
01593       InView = SYMBOLIC_VEC;
01594   if(List.isParameter("symminmultivec"))
01595     if(Teuchos::getParameter<bool>(List, "symminmultivec"))
01596       InView = SYMBOLIC_MULTIVEC;
01597   if(List.isParameter("symmoutvec"))
01598     if(Teuchos::getParameter<bool>(List, "symmoutvec"))
01599       OutView = SYMBOLIC_VEC;
01600   if(List.isParameter("symmoutmultivec"))
01601     if(Teuchos::getParameter<bool>(List, "symmoutmultivec"))
01602       OutView = SYMBOLIC_MULTIVEC;
01603   if(List.isParameter("symmintervec"))
01604     if(Teuchos::getParameter<bool>(List, "symmintermultivec"))
01605       IntermediateView = SYMBOLIC_MULTIVEC;
01606   if(List.isParameter("invalidinter"))
01607     if(Teuchos::getParameter<bool>(List, "invalidinter"))
01608       IntermediateView = NULL;
01609   if(this->Comm().NumProc() > 1) {
01610     if(InVec.NumVectors() == 1)
01611       InView = SYMBOLIC_VEC;
01612     else
01613       InView = SYMBOLIC_MULTIVEC;
01614     if(OutVec.NumVectors() == 1)
01615       OutView = SYMBOLIC_VEC;
01616     else
01617       OutView = SYMBOLIC_MULTIVEC;
01618     if(Intermediate.NumVectors() == 1)
01619       IntermediateView = SYMBOLIC_VEC;
01620     else
01621       IntermediateView = SYMBOLIC_MULTIVEC;
01622   }
01623   if(TransA)
01624     ReturnVal = oski_SetHintMatPowMult(A_tunable_, OP_TRANS, Power, Alpha, InView, Beta, OutView, IntermediateView, NumCalls);
01625   else
01626     ReturnVal = oski_SetHintMatPowMult(A_tunable_, OP_NORMAL, Power, Alpha, InView, Beta, OutView, IntermediateView, NumCalls);
01627   if(ReturnVal)
01628     std::cerr << "Set hint matpow multiply error\n";
01629   return ReturnVal;
01630 }
01631 
01632 
01633 int Epetra_OskiMatrix::TuneMatrix() {
01634   int ReturnVal;
01635   ReturnVal = oski_TuneMat(A_tunable_);
01636   if(IsMatrixTransformed())
01637     Copy_Created_ = true;
01638   return ReturnVal;
01639 }
01640 
01641 int Epetra_OskiMatrix::IsMatrixTransformed() const {
01642   return oski_IsMatPermuted(A_tunable_);
01643 }
01644 
01645 const Epetra_OskiMatrix& Epetra_OskiMatrix::ViewTransformedMat() const {
01646   //might need a more efficient way to do this
01647   Epetra_OskiMatrix* Transformed = NULL;
01648   Epetra_OskiMatrix Temp(*this); //should be some lightweight copy
01649   Transformed = &Temp;
01650   Transformed->A_tunable_ = const_cast <oski_matrix_t> (oski_ViewPermutedMat(A_tunable_));
01651   return *Transformed;
01652 }
01653 
01654 const Epetra_OskiPermutation& Epetra_OskiMatrix::ViewRowPermutation() const {
01655   Epetra_OskiPermutation* RowPerm = NULL;
01656   RowPerm = new Epetra_OskiPermutation[1];
01657   (*RowPerm).ReplacePermutation(const_cast <oski_perm_t> (oski_ViewPermutedMatRowPerm(A_tunable_)));
01658   return *RowPerm;
01659 }
01660 
01661 const Epetra_OskiPermutation& Epetra_OskiMatrix::ViewColumnPermutation() const {
01662   Epetra_OskiPermutation* ColPerm;
01663   ColPerm = new Epetra_OskiPermutation[1];
01664   (*ColPerm).ReplacePermutation(const_cast <oski_perm_t> (oski_ViewPermutedMatColPerm(A_tunable_)));
01665   return *ColPerm;
01666 }
01667 
01668 char* Epetra_OskiMatrix::GetMatrixTransforms() const {
01669   char* ReturnVal = NULL;
01670   ReturnVal = oski_GetMatTransforms(A_tunable_);
01671   if(ReturnVal == NULL)
01672     std::cerr << "Error in GetMatrixTransforms\n";
01673   return ReturnVal;
01674 }
01675 
01676 int Epetra_OskiMatrix::ApplyMatrixTransforms(const char* Transforms) {
01677   int ReturnVal;
01678   ReturnVal = oski_ApplyMatTransforms(A_tunable_, Transforms);
01679   if(ReturnVal)
01680     std::cerr << "Error in ApplyMatrixTransforms\n";
01681   return ReturnVal;
01682 }
01683 #endif
01684 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines