Epetra_OskiMatrix.cpp

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

Generated on Wed May 12 21:41:05 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7