EpetraExt Package Browser (Single Doxygen Collection) Development
hypre_Helpers.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include "hypre_Helpers.hpp"
00043 #include "Teuchos_UnitTestHarness.hpp"
00044 #include "Teuchos_Assert.hpp"
00045 #include "Teuchos_Array.hpp"
00046 #include <iostream>
00047 #include <fstream>
00048 
00049 #include <math.h>
00050 #include <cstdlib>
00051 #include <time.h>
00052 #include <vector>
00053 #include <map>
00054 
00055 EpetraExt_HypreIJMatrix::EpetraExt_HypreIJMatrix* newHypreMatrix(const int N, const int type)
00056 {
00057   HYPRE_IJMatrix Matrix;
00058   int ierr = 0;
00059   int i;
00060   int numprocs, rank;
00061   MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
00062   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00063   
00064   int ilower = (N/numprocs)*rank;
00065   int iupper = (N/numprocs)*(rank+1);
00066   if(rank == numprocs-1){ /*Last processor */
00067     iupper = N;
00068   }
00069   
00070   // Create
00071   ierr += HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper-1, ilower, iupper-1, &Matrix);
00072   ierr += HYPRE_IJMatrixSetObjectType(Matrix,HYPRE_PARCSR);
00073   // Initialize
00074   ierr += HYPRE_IJMatrixInitialize(Matrix);
00075   
00076   if(Matrix == NULL){
00077     printf("Error! Matrix is NULL!\n");
00078     std::cin >> ierr;
00079   }
00080   
00081   if(type == 0){
00082     // Set values
00083     int rows[1];
00084     int cols[1];
00085     double values[1];
00086     int ncols = 1;
00087     for(i = ilower; i < iupper; i++) {
00088       rows[0] = i;
00089       cols[0] = i;
00090       values[0] = 1.0;
00091       ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, cols, values);
00092     }
00093   } else if(type == 1){
00094     srand(time(NULL));
00095     // Set values
00096     int rows[1];
00097     int cols[1];
00098     double values[1];
00099     int ncols = 1;
00100     for(i = ilower; i < iupper; i++) {
00101       rows[0] = i;
00102       cols[0] = i;
00103       values[0] =  ( (double)rand()/(double)RAND_MAX ) * 100;
00104       ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, cols, values);
00105     }
00106     
00107   } else if(type == 2){
00108     // Set values
00109     int rows[1];
00110     Teuchos::Array<int> cols; cols.resize(N);
00111     Teuchos::Array<double> values; values.resize(N);
00112     int ncols = N;
00113     for(i = ilower; i < iupper; i++) {
00114       for(int j = 0; j < N; j++){
00115         cols[j] = j;
00116         values[j] = j;
00117       }
00118       rows[0] = i;
00119       ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
00120     }
00121   } else if(type == 3){
00122     srand(time(NULL));
00123     int rows[1];
00124     Teuchos::Array<int> cols; cols.resize(N);
00125     Teuchos::Array<double> values; values.resize(N);
00126     int ncols = N;
00127     for(i = ilower; i < iupper; i++) {
00128       for(int j = 0; j < N; j++){
00129         cols[j] = j;
00130         double currVal =  ( (double)rand()/(double)RAND_MAX ) * 100;
00131         values[j] = currVal;
00132       }
00133       rows[0] = i;
00134       ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
00135     }
00136   } else {
00137     srand(time(NULL));
00138     int rows[1];
00139     for(i = ilower; i < iupper; i++) {
00140       int ncols = (int)(1+( (double)rand()/(double)RAND_MAX ) * 0.5*(N-1));
00141       TEUCHOS_TEST_FOR_EXCEPTION(ncols <= 0, std::logic_error, "ncols is negative");
00142       Teuchos::Array<int> cols; cols.resize(ncols);
00143       Teuchos::Array<double> values; values.resize(ncols);
00144       for(int j = 0; j < ncols; j++){
00145         int index = 0;
00146         if(i-(ncols/2) >= 0 && i+(ncols/2) < N){
00147           index = j + i - (ncols/2);
00148         } else if (i-(ncols/2) < 0){
00149           index = j + i;
00150         } else{
00151           index = j + i - (ncols-1);
00152           }
00153         
00154         cols[j] = index;
00155         double currVal =  ( (double)rand()/(double)RAND_MAX ) * 100;
00156         values[j] = currVal;
00157       }
00158       rows[0] = i;
00159       ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
00160     }
00161   }
00162   // Assemble
00163   ierr += HYPRE_IJMatrixAssemble(Matrix);
00164   EpetraExt_HypreIJMatrix* RetMat = new EpetraExt_HypreIJMatrix(Matrix);
00165   //Don't HYPRE_IJMatrixDestroy(Matrix);
00166   return RetMat;
00167 }
00168 
00169 Epetra_CrsMatrix::Epetra_CrsMatrix* newCrsMatrix(EpetraExt_HypreIJMatrix &Matrix)
00170 {
00171   int N = Matrix.NumGlobalRows();
00172   Epetra_CrsMatrix* TestMat = new Epetra_CrsMatrix(Copy, Matrix.RowMatrixRowMap(), Matrix.RowMatrixColMap(), N, false);
00173   
00174   for(int i = 0; i < Matrix.NumMyRows(); i++){
00175     int entries;
00176     Matrix.NumMyRowEntries(i,entries);
00177     Teuchos::Array<double> Values; Values.resize(entries);
00178     Teuchos::Array<int> Indices; Indices.resize(entries);
00179     int NumEntries;
00180     Matrix.ExtractMyRowCopy(i,entries,NumEntries,&Values[0], &Indices[0]);
00181     for(int j = 0; j < NumEntries; j++){
00182       double currVal[1];
00183       currVal[0] = Values[j];
00184       int currInd[1];
00185       currInd[0] = Matrix.RowMatrixColMap().GID(Indices[j]);
00186       TestMat->InsertGlobalValues(Matrix.RowMatrixRowMap().GID(i), 1, currVal, currInd);
00187     }
00188   }
00189   TestMat->FillComplete();
00190   return TestMat;
00191 }
00192 
00193 bool EquivalentVectors(Epetra_MultiVector &Y1, Epetra_MultiVector &Y2, const double tol){
00194   
00195   bool retVal = true;
00196   
00197   int num_vectors = Y1.NumVectors();
00198   if(Y2.NumVectors() != num_vectors){
00199     printf("Multivectors do not have same number of vectors.\n");
00200     return false;
00201   }
00202   
00203   for(int j = 0; j < num_vectors; j++){
00204     if(Y1.MyLength() != Y2.MyLength()){
00205       printf("Vectors are not same size on local processor.\n");
00206       return false;
00207     }
00208     Teuchos::Array<double> Y1_vals; Y1_vals.resize(Y1.MyLength());
00209     Teuchos::Array<double> Y2_vals; Y2_vals.resize(Y2.MyLength());
00210     (*Y1(j)).ExtractCopy(&Y1_vals[0]);
00211     (*Y2(j)).ExtractCopy(&Y2_vals[0]);
00212     
00213     for(int i = 0; i < Y1.MyLength(); i++){
00214       if(fabs(Y1_vals[i] - Y2_vals[i]) > tol){
00215         printf("Vector number[%d] ", j);
00216         printf("Val1[%d] = %f != Val2[%d] = %f\n", i, Y1_vals[i], i, Y2_vals[i]);
00217         retVal = false;
00218       }
00219     }
00220   }
00221   if(retVal == false){
00222     printf("Failed equivalent vectors.\n");
00223   }
00224   return retVal;
00225 }
00226 
00227 
00228 
00229 bool EquivalentMatrices(Epetra_RowMatrix &HypreMatrix, Epetra_RowMatrix &CrsMatrix, const double tol){
00230   bool retVal = true;
00231   for(int j = 0; j < HypreMatrix.NumMyRows(); j++){
00232     
00233     int NumEntries;
00234     int entries1;
00235     int entries2;
00236     HypreMatrix.NumMyRowEntries(j,NumEntries);
00237 
00238     Teuchos::Array<double> Y1_vals; Y1_vals.resize(NumEntries);
00239     Teuchos::Array<double> Y2_vals; Y2_vals.resize(NumEntries);
00240     Teuchos::Array<int> indices1; indices1.resize(NumEntries);
00241     Teuchos::Array<int> indices2; indices2.resize(NumEntries);
00242      
00243     HypreMatrix.ExtractMyRowCopy(j,NumEntries, entries1, &Y1_vals[0], &indices1[0]);
00244     CrsMatrix.ExtractMyRowCopy(j,NumEntries, entries2, &Y2_vals[0], &indices2[0]);
00245 
00246     std::map<int,double> Y1map;
00247     std::map<int,double> Y2map;
00248     for (int i=0; i < NumEntries ; ++i) {
00249       Y1map[indices1[i]] = Y1_vals[i]; 
00250       Y2map[indices2[i]] = Y2_vals[i]; 
00251     }
00252     retVal = retVal && (Y1map == Y2map);
00253   }
00254   Teuchos::Array<int> vals; vals.resize(HypreMatrix.Comm().NumProc());
00255   int my_vals[1]; my_vals[0] = (int)retVal;
00256   HypreMatrix.Comm().GatherAll(my_vals, &vals[0], 1);
00257   for(int i = 0; i < HypreMatrix.Comm().NumProc(); i++){
00258     if(vals[i] == false){
00259       retVal = false;
00260     }
00261   }
00262   if(retVal == false){
00263     printf("[%d]Failed matrix equivalency test.\n", HypreMatrix.Comm().MyPID());
00264     if(HypreMatrix.Comm().MyPID() == 0){
00265       //std::ofstream outfile("Matrices.txt");
00266       
00267       //HypreMatrix.Print(outfile);
00268       //CrsMatrix.Print(outfile);
00269       //outfile.close();
00270     }
00271   }
00272   return retVal;
00273 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines