rectMatrixTest.cc

Go to the documentation of this file.
00001 //@HEADER
00002 /*
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 */
00029 //@HEADER
00030 
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include <ctype.h>
00034 #include <assert.h>
00035 #include <string.h>
00036 #include <math.h>
00037 #include "Petra_Comm.h"
00038 #include "Petra_Map.h"
00039 #include "Petra_Time.h"
00040 #include "Petra_RDP_MultiVector.h"
00041 #include "Petra_RDP_Vector.h"
00042 #include "Petra_RDP_CRS_Matrix.h"
00043 #ifdef PETRA_MPI
00044 #include "mpi.h"
00045 #endif
00046 #ifndef __cplusplus
00047 #define __cplusplus
00048 #endif
00049 
00050 // Test code to make a rectangular petra matrix from data in a file
00051 // -- vh
00052 
00053 // prototypes
00054 Petra_RDP_CRS_Matrix* readMatrixIn(FILE *dataFile, Petra_Comm& Comm);
00055 Petra_RDP_CRS_Matrix* readRectMatrixIn(FILE *dataFile, Petra_Comm& Comm);
00056 Petra_RDP_Vector* readVectorIn(FILE *dataFile, Petra_Comm& Comm);
00057 void matVecTest(Petra_RDP_CRS_Matrix* Aptr,
00058                 Petra_RDP_Vector* xptr,
00059                 Petra_RDP_Vector* bptr);
00060 
00061 
00062 int main(int argc, char *argv[])
00063 {
00064   int ierr = 0, i, j;
00065   
00066 #ifdef PETRA_MPI
00067   // Initialize MPI
00068   MPI_Init(&argc,&argv);
00069   int size, rank; // Number of MPI processes, My process ID
00070   MPI_Comm_size(MPI_COMM_WORLD, &size);
00071   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00072 #else
00073   int size = 1; // Serial case (not using MPI)
00074   int rank = 0;
00075 #endif
00076 
00077 #ifdef PETRA_MPI
00078   Petra_Comm & Comm = *new Petra_Comm( MPI_COMM_WORLD );
00079 #else
00080   Petra_Comm & Comm = *new Petra_Comm();
00081 #endif
00082 
00083   int MyPID = Comm.MyPID();
00084   int NumProc = Comm.NumProc();
00085   // cout << "Processor "<<MyPID<<" of "<< NumProc << " is alive."<<endl;
00086 
00087   bool verbose = (MyPID==0);
00088 
00089 
00090   // read in matrices:
00091   FILE *dataFile;
00092 
00093   // This is a square matrix
00094   cout << " reading in F matrix " << endl;
00095   dataFile = fopen("F.data","r");
00096   Petra_RDP_CRS_Matrix* Fptr = readMatrixIn(dataFile, Comm);
00097   fclose(dataFile);
00098 
00099   // Read in my rectangular matrix
00100   cout << " reading in B matrix " << endl;
00101   dataFile = fopen("B.data","r");
00102   Petra_RDP_CRS_Matrix* Bptr = readRectMatrixIn(dataFile, Comm);
00103   fclose(dataFile);
00104   Petra_RDP_CRS_Matrix& B = *Bptr;
00105   cout << "global rows " << B.NumGlobalRows() << endl;
00106   cout << "global cols " << B.NumGlobalCols() << endl;
00107   cout << "my local rows " << B.NumMyRows() << endl;
00108   cout << "my local cols  " << B.NumMyCols() << endl;
00109 
00110   // read in vectors (b and soln)
00111   cout << "reading in vector b " << endl;
00112   dataFile = fopen("rhs.data","r");
00113   Petra_RDP_Vector* bptr = readVectorIn(dataFile, Comm);
00114   fclose(dataFile);
00115 
00116 }
00117 
00118 
00119 Petra_RDP_CRS_Matrix* readMatrixIn(FILE *dataFile, Petra_Comm& Comm)
00120 {
00128 
00129   int NumGlobalElements = 0;
00130   int maxNumNz = 0;
00131   int i, j;
00132   // dataFile = fopen("B.data","r");
00133   fscanf(dataFile, "%d", &NumGlobalElements);
00134   // cout << "NumGlobalElements = " << NumGlobalElements << "\n" << endl;
00135   fscanf(dataFile, "%d", &maxNumNz);
00136 
00137   // Construct a Map that puts approximately the same number of 
00138   // equations on each processor.
00139   Petra_Map& Map = *new Petra_Map(NumGlobalElements, 0, Comm);
00140 
00141   // Get update list and number of local equations from newly created Map.
00142   int NumMyElements = Map.NumMyElements();
00143 
00144   int * MyGlobalElements = new int[NumMyElements];
00145   Map.MyGlobalElements(MyGlobalElements);
00146 
00147   int * NumNz = new int[NumMyElements];
00148   for (i=0; i<NumMyElements; i++)
00149     {
00150       fscanf(dataFile, "%d", &NumNz[i]);
00151       // NumNz[i] = NumNz[i] - 1; // subtracting off one for the diagonals
00152       // figure out what to do for this in the non-square matrices (B)
00153       // cout << NumNz[i] << endl;
00154     }
00155 
00156   Petra_RDP_CRS_Matrix* Mptr = new Petra_RDP_CRS_Matrix(Copy, Map, NumNz);
00157   Petra_RDP_CRS_Matrix& M = *Mptr;
00158   double *Values = new double[maxNumNz];
00159   int *Indices = new int[maxNumNz];
00160   int NumEntries;
00164   int nnzM = 0; 
00165   
00166   fscanf(dataFile, "%d", &nnzM);
00167   // cout << "\nnumber of nonzeros in B: " << nnzB << "\n" << endl;
00168 
00169 
00170   int * iM = new int[nnzM];
00171   int * jM = new int[nnzM];
00172   double * sM = new double[nnzM];
00173   for (i=0; i<nnzM; i++)
00174     {
00175       fscanf(dataFile, "%d %d %lg", &iM[i], &jM[i], &sM[i]);
00176       // matlab used indexing from 1, C++ starts at 0
00177       // so subtract 1:
00178       iM[i]--;
00179       jM[i]--;
00180       // cout << "iM: " << iM[i]
00181       // << "\tjM: " << jM[i]
00182       // << "\tsM: " << sM[i]
00183       // << endl;
00184     }
00185 
00187   int offset = 0;
00188   for (i=0; i<NumGlobalElements; i++) 
00189     {
00190       // cout << "NumNz[" << i << "] = " << NumNz[i] << endl;
00191       for (j=0; j<NumNz[i]; j++)
00192         {
00194           Indices[j] = jM[offset + j];
00195           Values[j] = sM[offset + j];
00196           // cout << "iM: " << iM[offset + j]
00197           //     << "\tjM: " << jM[offset + j]
00198           //     << "\tsM: " << sM[offset + j]
00199           //     << endl;
00200         }
00201       NumEntries = NumNz[i];
00202       assert(M.InsertGlobalValues(MyGlobalElements[i], 
00203                 NumEntries, Values, Indices)==0);
00204       // cout << "offset = " << offset << endl;
00205       offset = offset + NumNz[i];
00206       // cout << "Got to here in B matrix." <<endl;
00207     }
00208 
00209   // Finish up
00210   assert(M.FillComplete()==0);
00211 
00212   cout << "nonzeros = " << M.NumGlobalNonzeros() << endl;
00213   cout << "rows = " << M.NumGlobalRows() << endl;
00214   cout << "cols = " << M.NumGlobalCols() << endl;
00215   return Mptr;
00216 }
00217 
00218 
00219 Petra_RDP_CRS_Matrix* readRectMatrixIn(FILE *dataFile, Petra_Comm& Comm)
00220 {
00228 
00229   int NumGlobalElements = 0;
00230   int maxNumNz = 0;
00231   int i, j;
00232   // dataFile = fopen("B.data","r");
00233   fscanf(dataFile, "%d", &NumGlobalElements);
00234   // cout << "NumGlobalElements = " << NumGlobalElements << "\n" << endl;
00235   fscanf(dataFile, "%d", &maxNumNz);
00236 
00237   // Construct a Map that puts approximately the same number of 
00238   // equations on each processor.
00239   Petra_Map& Map = *new Petra_Map(NumGlobalElements, 0, Comm);
00240 
00241   // make another map for the columns'
00242   // Note -- figure out how to pass in the number of columns
00243   Petra_Map& ColMap = *new Petra_Map(24, 0, Comm);
00244 
00245   // Get update list and number of local equations from newly created Map.
00246   int NumMyElements = Map.NumMyElements();
00247 
00248   int * MyGlobalElements = new int[NumMyElements];
00249   Map.MyGlobalElements(MyGlobalElements);
00250 
00251   int * NumNz = new int[NumMyElements];
00252   for (i=0; i<NumMyElements; i++)
00253     {
00254       fscanf(dataFile, "%d", &NumNz[i]);
00255       // NumNz[i] = NumNz[i] - 1; // subtracting off one for the diagonals
00256       // figure out what to do for this in the non-square matrices (B)
00257       // cout << NumNz[i] << endl;
00258     }
00259 
00260   Petra_RDP_CRS_Matrix* Mptr = new Petra_RDP_CRS_Matrix(Copy, Map, NumNz);
00261   Petra_RDP_CRS_Matrix& M = *Mptr;
00262   double *Values = new double[maxNumNz];
00263   int *Indices = new int[maxNumNz];
00264   int NumEntries;
00268   int nnzM = 0; 
00269   
00270   fscanf(dataFile, "%d", &nnzM);
00271   // cout << "\nnumber of nonzeros in B: " << nnzB << "\n" << endl;
00272 
00273 
00274   int * iM = new int[nnzM];
00275   int * jM = new int[nnzM];
00276   double * sM = new double[nnzM];
00277   for (i=0; i<nnzM; i++)
00278     {
00279       fscanf(dataFile, "%d %d %lg", &iM[i], &jM[i], &sM[i]);
00280       // matlab used indexing from 1, C++ starts at 0
00281       // so subtract 1:
00282       iM[i]--;
00283       jM[i]--;
00284       // cout << "iM: " << iM[i]
00285       // << "\tjM: " << jM[i]
00286       // << "\tsM: " << sM[i]
00287       // << endl;
00288     }
00289 
00291   int offset = 0;
00292   for (i=0; i<NumGlobalElements; i++) 
00293     {
00294       // cout << "NumNz[" << i << "] = " << NumNz[i] << endl;
00295       for (j=0; j<NumNz[i]; j++)
00296         {
00298           Indices[j] = jM[offset + j];
00299           Values[j] = sM[offset + j];
00300           // cout << "iM: " << iM[offset + j]
00301           //     << "\tjM: " << jM[offset + j]
00302           //     << "\tsM: " << sM[offset + j]
00303           //     << endl;
00304         }
00305       NumEntries = NumNz[i];
00306       assert(M.InsertGlobalValues(MyGlobalElements[i], 
00307                 NumEntries, Values, Indices)==0);
00308       // cout << "offset = " << offset << endl;
00309       offset = offset + NumNz[i];
00310       // cout << "Got to here in B matrix." <<endl;
00311     }
00312 
00313   // Finish up
00314   assert(M.FillComplete(ColMap, Map)==0);
00315 
00316   cout << "nonzeros = " << M.NumGlobalNonzeros() << endl;
00317   cout << "rows = " << M.NumGlobalRows() << endl;
00318   cout << "cols = " << M.NumGlobalCols() << endl;
00319   return Mptr;
00320 }
00321 
00322 
00323 
00324 
00325 
00326 Petra_RDP_Vector* readVectorIn(FILE *dataFile, Petra_Comm& Comm)
00327 {
00333 
00334   int NumGlobalElements = 0;
00335   int NzElms = 0;
00336   int i;
00337 
00338   fscanf(dataFile, "%d", &NumGlobalElements);
00339   fscanf(dataFile, "%d", &NzElms);
00340   // Construct a Map that puts approximately the same number of 
00341   // equations on each processor.
00342   Petra_Map& Map = *new Petra_Map(NumGlobalElements, 0, Comm);
00343 
00344   // Get update list and number of local equations from newly created Map.
00345   int NumMyElements = Map.NumMyElements();
00346   int * MyGlobalElements = new int[NumMyElements];
00347   Map.MyGlobalElements(MyGlobalElements);
00348 
00349   // make a petra map filled with zeros
00350   Petra_RDP_Vector* vptr = new Petra_RDP_Vector(Map);
00351   Petra_RDP_Vector& v = *vptr;
00352 
00353   // now fill in the nonzero elements
00354   double * myArray = new double[NumMyElements];
00355   int tempInd;
00356   double tempVal;
00357   // cout << "Length v " << NumGlobalElements << endl;
00358   // cout << "NzElms " << NzElms << endl;
00359   for (i=0; i<NzElms; i++)
00360     {
00361       fscanf(dataFile, "%d %lg", &tempInd, &tempVal);
00362       v[tempInd] = tempVal;
00363       // cout << tempVal << endl;
00364     }
00365   //  Petra_RDP_CRS_Matrix& M = *new Petra_RDP_CRS_Matrix(Copy, Map, NumNz);
00366  
00367   return vptr;
00368 
00369 }
00370 
00371 // small matrix vector multiply test
00372 void matVecTest(Petra_RDP_CRS_Matrix* Aptr,
00373                 Petra_RDP_Vector* xptr,
00374                 Petra_RDP_Vector* bptr)
00375 {
00376   Petra_RDP_CRS_Matrix& A = *Aptr;
00377   Petra_RDP_Vector& x = *xptr;
00378   Petra_RDP_Vector& b = *bptr; // we're going to overwrite b anyway
00379   // will that overwrite x, too? Look at ExtractCopy for alternative.
00380 
00381   A.Multiply(1, x, b);
00382 }
00383 
00384 
00385 
00386 // matrix vector multiply for our block matrix
00387 /************
00388 Petra_RDP_Vector* = myMatVecMult(Petra_RDP_CRS_Matrix* Bptr, 
00389                                  Petra_RDP_CRS_Matrix* Fptr,
00390                                  Petra_RDP_CRS_Matrix* Cptr,
00391                                  Petra_RDP_Vector* xptr)
00392 
00393 {
00394   // A = [F B'; B C]
00395   // return Ax pointer
00396   // cout << "block matrix-vector multiply" << endl;
00397 
00398 }
00399 *************/

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1