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