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 (2009) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #include "hypre_Helpers.hpp"
00030 #include "Teuchos_UnitTestHarness.hpp"
00031 #include "Teuchos_TestForException.hpp"
00032 #include "Teuchos_Array.hpp"
00033 #include <iostream>
00034 #include <fstream>
00035 
00036 #include <math.h>
00037 #include <stdlib.h>
00038 #include <time.h>
00039 #include <vector>
00040 #include <map>
00041 
00042 EpetraExt_HypreIJMatrix::EpetraExt_HypreIJMatrix* newHypreMatrix(const int N, const int type)
00043 {
00044   HYPRE_IJMatrix Matrix;
00045   int ierr = 0;
00046   int i;
00047   int numprocs, rank;
00048   MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
00049   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00050   
00051   int ilower = (N/numprocs)*rank;
00052   int iupper = (N/numprocs)*(rank+1);
00053   if(rank == numprocs-1){ /*Last processor */
00054     iupper = N;
00055   }
00056   
00057   // Create
00058   ierr += HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper-1, ilower, iupper-1, &Matrix);
00059   ierr += HYPRE_IJMatrixSetObjectType(Matrix,HYPRE_PARCSR);
00060   // Initialize
00061   ierr += HYPRE_IJMatrixInitialize(Matrix);
00062   
00063   if(Matrix == NULL){
00064     printf("Error! Matrix is NULL!\n");
00065     std::cin >> ierr;
00066   }
00067   
00068   if(type == 0){
00069     // Set values
00070     int rows[1];
00071     int cols[1];
00072     double values[1];
00073     int ncols = 1;
00074     for(i = ilower; i < iupper; i++) {
00075       rows[0] = i;
00076       cols[0] = i;
00077       values[0] = 1.0;
00078       ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, cols, values);
00079     }
00080   } else if(type == 1){
00081     srand(time(NULL));
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] =  ( (double)rand()/(double)RAND_MAX ) * 100;
00091       ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, cols, values);
00092     }
00093     
00094   } else if(type == 2){
00095     // Set values
00096     int rows[1];
00097     Teuchos::Array<int> cols; cols.resize(N);
00098     Teuchos::Array<double> values; values.resize(N);
00099     int ncols = N;
00100     for(i = ilower; i < iupper; i++) {
00101       for(int j = 0; j < N; j++){
00102         cols[j] = j;
00103         values[j] = j;
00104       }
00105       rows[0] = i;
00106       ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
00107     }
00108   } else if(type == 3){
00109     srand(time(NULL));
00110     int rows[1];
00111     Teuchos::Array<int> cols; cols.resize(N);
00112     Teuchos::Array<double> values; values.resize(N);
00113     int ncols = N;
00114     for(i = ilower; i < iupper; i++) {
00115       for(int j = 0; j < N; j++){
00116         cols[j] = j;
00117         double currVal =  ( (double)rand()/(double)RAND_MAX ) * 100;
00118         values[j] = currVal;
00119       }
00120       rows[0] = i;
00121       ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
00122     }
00123   } else {
00124     srand(time(NULL));
00125     int rows[1];
00126     for(i = ilower; i < iupper; i++) {
00127       int ncols = (int)(1+( (double)rand()/(double)RAND_MAX ) * 0.5*(N-1));
00128       TEST_FOR_EXCEPTION(ncols <= 0, std::logic_error, "ncols is negative");
00129       Teuchos::Array<int> cols; cols.resize(ncols);
00130       Teuchos::Array<double> values; values.resize(ncols);
00131       for(int j = 0; j < ncols; j++){
00132         int index = 0;
00133         if(i-(ncols/2) >= 0 && i+(ncols/2) < N){
00134           index = j + i - (ncols/2);
00135         } else if (i-(ncols/2) < 0){
00136           index = j + i;
00137         } else{
00138           index = j + i - (ncols-1);
00139           }
00140         
00141         cols[j] = index;
00142         double currVal =  ( (double)rand()/(double)RAND_MAX ) * 100;
00143         values[j] = currVal;
00144       }
00145       rows[0] = i;
00146       ierr += HYPRE_IJMatrixSetValues(Matrix, 1, &ncols, rows, &cols[0], &values[0]);
00147     }
00148   }
00149   // Assemble
00150   ierr += HYPRE_IJMatrixAssemble(Matrix);
00151   EpetraExt_HypreIJMatrix* RetMat = new EpetraExt_HypreIJMatrix(Matrix);
00152   //Don't HYPRE_IJMatrixDestroy(Matrix);
00153   return RetMat;
00154 }
00155 
00156 Epetra_CrsMatrix::Epetra_CrsMatrix* newCrsMatrix(EpetraExt_HypreIJMatrix &Matrix)
00157 {
00158   int N = Matrix.NumGlobalRows();
00159   Epetra_CrsMatrix* TestMat = new Epetra_CrsMatrix(Copy, Matrix.RowMatrixRowMap(), Matrix.RowMatrixColMap(), N, false);
00160   
00161   for(int i = 0; i < Matrix.NumMyRows(); i++){
00162     int entries;
00163     Matrix.NumMyRowEntries(i,entries);
00164     Teuchos::Array<double> Values; Values.resize(entries);
00165     Teuchos::Array<int> Indices; Indices.resize(entries);
00166     int NumEntries;
00167     Matrix.ExtractMyRowCopy(i,entries,NumEntries,&Values[0], &Indices[0]);
00168     for(int j = 0; j < NumEntries; j++){
00169       double currVal[1];
00170       currVal[0] = Values[j];
00171       int currInd[1];
00172       currInd[0] = Matrix.RowMatrixColMap().GID(Indices[j]);
00173       TestMat->InsertGlobalValues(Matrix.RowMatrixRowMap().GID(i), 1, currVal, currInd);
00174     }
00175   }
00176   TestMat->FillComplete();
00177   return TestMat;
00178 }
00179 
00180 bool EquivalentVectors(Epetra_MultiVector &Y1, Epetra_MultiVector &Y2, const double tol){
00181   
00182   bool retVal = true;
00183   
00184   int num_vectors = Y1.NumVectors();
00185   if(Y2.NumVectors() != num_vectors){
00186     printf("Multivectors do not have same number of vectors.\n");
00187     return false;
00188   }
00189   
00190   for(int j = 0; j < num_vectors; j++){
00191     if(Y1.MyLength() != Y2.MyLength()){
00192       printf("Vectors are not same size on local processor.\n");
00193       return false;
00194     }
00195     Teuchos::Array<double> Y1_vals; Y1_vals.resize(Y1.MyLength());
00196     Teuchos::Array<double> Y2_vals; Y2_vals.resize(Y2.MyLength());
00197     (*Y1(j)).ExtractCopy(&Y1_vals[0]);
00198     (*Y2(j)).ExtractCopy(&Y2_vals[0]);
00199     
00200     for(int i = 0; i < Y1.MyLength(); i++){
00201       if(fabs(Y1_vals[i] - Y2_vals[i]) > tol){
00202         printf("Vector number[%d] ", j);
00203         printf("Val1[%d] = %f != Val2[%d] = %f\n", i, Y1_vals[i], i, Y2_vals[i]);
00204         retVal = false;
00205       }
00206     }
00207   }
00208   if(retVal == false){
00209     printf("Failed equivalent vectors.\n");
00210   }
00211   return retVal;
00212 }
00213 
00214 
00215 
00216 bool EquivalentMatrices(Epetra_RowMatrix &HypreMatrix, Epetra_RowMatrix &CrsMatrix, const double tol){
00217   bool retVal = true;
00218   for(int j = 0; j < HypreMatrix.NumMyRows(); j++){
00219     
00220     int NumEntries;
00221     int entries1;
00222     int entries2;
00223     HypreMatrix.NumMyRowEntries(j,NumEntries);
00224 
00225     Teuchos::Array<double> Y1_vals; Y1_vals.resize(NumEntries);
00226     Teuchos::Array<double> Y2_vals; Y2_vals.resize(NumEntries);
00227     Teuchos::Array<int> indices1; indices1.resize(NumEntries);
00228     Teuchos::Array<int> indices2; indices2.resize(NumEntries);
00229      
00230     HypreMatrix.ExtractMyRowCopy(j,NumEntries, entries1, &Y1_vals[0], &indices1[0]);
00231     CrsMatrix.ExtractMyRowCopy(j,NumEntries, entries2, &Y2_vals[0], &indices2[0]);
00232 
00233     std::map<int,double> Y1map;
00234     std::map<int,double> Y2map;
00235     for (int i=0; i < NumEntries ; ++i) {
00236       Y1map[indices1[i]] = Y1_vals[i]; 
00237       Y2map[indices2[i]] = Y2_vals[i]; 
00238     }
00239     retVal = retVal && (Y1map == Y2map);
00240   }
00241   Teuchos::Array<int> vals; vals.resize(HypreMatrix.Comm().NumProc());
00242   int my_vals[1]; my_vals[0] = (int)retVal;
00243   HypreMatrix.Comm().GatherAll(my_vals, &vals[0], 1);
00244   for(int i = 0; i < HypreMatrix.Comm().NumProc(); i++){
00245     if(vals[i] == false){
00246       retVal = false;
00247     }
00248   }
00249   if(retVal == false){
00250     printf("[%d]Failed matrix equivalency test.\n", HypreMatrix.Comm().MyPID());
00251     if(HypreMatrix.Comm().MyPID() == 0){
00252       //std::ofstream outfile("Matrices.txt");
00253       
00254       //HypreMatrix.Print(outfile);
00255       //CrsMatrix.Print(outfile);
00256       //outfile.close();
00257     }
00258   }
00259   return retVal;
00260 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines