Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
rectMatrixTest.cc
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Epetra: Linear Algebra Services Package
5// Copyright 2011 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40//@HEADER
41
42
43#include <stdio.h>
44#include <stdlib.h>
45#include <ctype.h>
46#include <assert.h>
47#include <string.h>
48#include <math.h>
49#include "Petra_Comm.h"
50#include "Petra_Map.h"
51#include "Petra_Time.h"
52#include "Petra_RDP_MultiVector.h"
53#include "Petra_RDP_Vector.h"
54#include "Petra_RDP_CRS_Matrix.h"
55#ifdef PETRA_MPI
56#include "mpi.h"
57#endif
58#ifndef __cplusplus
59#define __cplusplus
60#endif
61
62// Test code to make a rectangular petra matrix from data in a file
63// -- vh
64
65// prototypes
66Petra_RDP_CRS_Matrix* readMatrixIn(FILE *dataFile, Petra_Comm& Comm);
67Petra_RDP_CRS_Matrix* readRectMatrixIn(FILE *dataFile, Petra_Comm& Comm);
68Petra_RDP_Vector* readVectorIn(FILE *dataFile, Petra_Comm& Comm);
69void matVecTest(Petra_RDP_CRS_Matrix* Aptr,
70 Petra_RDP_Vector* xptr,
71 Petra_RDP_Vector* bptr);
72
73
74int main(int argc, char *argv[])
75{
76 int ierr = 0, i, j;
77
78#ifdef PETRA_MPI
79 // Initialize MPI
80 MPI_Init(&argc,&argv);
81 int size, rank; // Number of MPI processes, My process ID
82 MPI_Comm_size(MPI_COMM_WORLD, &size);
83 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
84#else
85 int size = 1; // Serial case (not using MPI)
86 int rank = 0;
87#endif
88
89#ifdef PETRA_MPI
90 Petra_Comm & Comm = *new Petra_Comm( MPI_COMM_WORLD );
91#else
92 Petra_Comm & Comm = *new Petra_Comm();
93#endif
94
95 int MyPID = Comm.MyPID();
96 int NumProc = Comm.NumProc();
97 // cout << "Processor "<<MyPID<<" of "<< NumProc << " is alive."<<endl;
98
99 bool verbose = (MyPID==0);
100
101
102 // read in matrices:
103 FILE *dataFile;
104
105 // This is a square matrix
106 cout << " reading in F matrix " << endl;
107 dataFile = fopen("F.data","r");
108 Petra_RDP_CRS_Matrix* Fptr = readMatrixIn(dataFile, Comm);
109 fclose(dataFile);
110
111 // Read in my rectangular matrix
112 cout << " reading in B matrix " << endl;
113 dataFile = fopen("B.data","r");
114 Petra_RDP_CRS_Matrix* Bptr = readRectMatrixIn(dataFile, Comm);
115 fclose(dataFile);
116 Petra_RDP_CRS_Matrix& B = *Bptr;
117 cout << "global rows " << B.NumGlobalRows() << endl;
118 cout << "global cols " << B.NumGlobalCols() << endl;
119 cout << "my local rows " << B.NumMyRows() << endl;
120 cout << "my local cols " << B.NumMyCols() << endl;
121
122 // read in vectors (b and soln)
123 cout << "reading in vector b " << endl;
124 dataFile = fopen("rhs.data","r");
125 Petra_RDP_Vector* bptr = readVectorIn(dataFile, Comm);
126 fclose(dataFile);
127
128}
129
130
131Petra_RDP_CRS_Matrix* readMatrixIn(FILE *dataFile, Petra_Comm& Comm)
132{
140
141 int NumGlobalElements = 0;
142 int maxNumNz = 0;
143 int i, j;
144 // dataFile = fopen("B.data","r");
145 fscanf(dataFile, "%d", &NumGlobalElements);
146 // cout << "NumGlobalElements = " << NumGlobalElements << "\n" << endl;
147 fscanf(dataFile, "%d", &maxNumNz);
148
149 // Construct a Map that puts approximately the same number of
150 // equations on each processor.
151 Petra_Map& Map = *new Petra_Map(NumGlobalElements, 0, Comm);
152
153 // Get update list and number of local equations from newly created Map.
154 int NumMyElements = Map.NumMyElements();
155
156 int * MyGlobalElements = new int[NumMyElements];
157 Map.MyGlobalElements(MyGlobalElements);
158
159 int * NumNz = new int[NumMyElements];
160 for (i=0; i<NumMyElements; i++)
161 {
162 fscanf(dataFile, "%d", &NumNz[i]);
163 // NumNz[i] = NumNz[i] - 1; // subtracting off one for the diagonals
164 // figure out what to do for this in the non-square matrices (B)
165 // cout << NumNz[i] << endl;
166 }
167
168 Petra_RDP_CRS_Matrix* Mptr = new Petra_RDP_CRS_Matrix(Copy, Map, NumNz);
169 Petra_RDP_CRS_Matrix& M = *Mptr;
170 double *Values = new double[maxNumNz];
171 int *Indices = new int[maxNumNz];
172 int NumEntries;
176 int nnzM = 0;
177
178 fscanf(dataFile, "%d", &nnzM);
179 // cout << "\nnumber of nonzeros in B: " << nnzB << "\n" << endl;
180
181
182 int * iM = new int[nnzM];
183 int * jM = new int[nnzM];
184 double * sM = new double[nnzM];
185 for (i=0; i<nnzM; i++)
186 {
187 fscanf(dataFile, "%d %d %lg", &iM[i], &jM[i], &sM[i]);
188 // matlab used indexing from 1, C++ starts at 0
189 // so subtract 1:
190 iM[i]--;
191 jM[i]--;
192 // cout << "iM: " << iM[i]
193 // << "\tjM: " << jM[i]
194 // << "\tsM: " << sM[i]
195 // << endl;
196 }
197
199 int offset = 0;
200 for (i=0; i<NumGlobalElements; i++)
201 {
202 // cout << "NumNz[" << i << "] = " << NumNz[i] << endl;
203 for (j=0; j<NumNz[i]; j++)
204 {
206 Indices[j] = jM[offset + j];
207 Values[j] = sM[offset + j];
208 // cout << "iM: " << iM[offset + j]
209 // << "\tjM: " << jM[offset + j]
210 // << "\tsM: " << sM[offset + j]
211 // << endl;
212 }
213 NumEntries = NumNz[i];
214 assert(M.InsertGlobalValues(MyGlobalElements[i],
215 NumEntries, Values, Indices)==0);
216 // cout << "offset = " << offset << endl;
217 offset = offset + NumNz[i];
218 // cout << "Got to here in B matrix." <<endl;
219 }
220
221 // Finish up
222 assert(M.FillComplete()==0);
223
224 cout << "nonzeros = " << M.NumGlobalNonzeros() << endl;
225 cout << "rows = " << M.NumGlobalRows() << endl;
226 cout << "cols = " << M.NumGlobalCols() << endl;
227 return Mptr;
228}
229
230
231Petra_RDP_CRS_Matrix* readRectMatrixIn(FILE *dataFile, Petra_Comm& Comm)
232{
240
241 int NumGlobalElements = 0;
242 int maxNumNz = 0;
243 int i, j;
244 // dataFile = fopen("B.data","r");
245 fscanf(dataFile, "%d", &NumGlobalElements);
246 // cout << "NumGlobalElements = " << NumGlobalElements << "\n" << endl;
247 fscanf(dataFile, "%d", &maxNumNz);
248
249 // Construct a Map that puts approximately the same number of
250 // equations on each processor.
251 Petra_Map& Map = *new Petra_Map(NumGlobalElements, 0, Comm);
252
253 // make another map for the columns'
254 // Note -- figure out how to pass in the number of columns
255 Petra_Map& ColMap = *new Petra_Map(24, 0, Comm);
256
257 // Get update list and number of local equations from newly created Map.
258 int NumMyElements = Map.NumMyElements();
259
260 int * MyGlobalElements = new int[NumMyElements];
261 Map.MyGlobalElements(MyGlobalElements);
262
263 int * NumNz = new int[NumMyElements];
264 for (i=0; i<NumMyElements; i++)
265 {
266 fscanf(dataFile, "%d", &NumNz[i]);
267 // NumNz[i] = NumNz[i] - 1; // subtracting off one for the diagonals
268 // figure out what to do for this in the non-square matrices (B)
269 // cout << NumNz[i] << endl;
270 }
271
272 Petra_RDP_CRS_Matrix* Mptr = new Petra_RDP_CRS_Matrix(Copy, Map, NumNz);
273 Petra_RDP_CRS_Matrix& M = *Mptr;
274 double *Values = new double[maxNumNz];
275 int *Indices = new int[maxNumNz];
276 int NumEntries;
280 int nnzM = 0;
281
282 fscanf(dataFile, "%d", &nnzM);
283 // cout << "\nnumber of nonzeros in B: " << nnzB << "\n" << endl;
284
285
286 int * iM = new int[nnzM];
287 int * jM = new int[nnzM];
288 double * sM = new double[nnzM];
289 for (i=0; i<nnzM; i++)
290 {
291 fscanf(dataFile, "%d %d %lg", &iM[i], &jM[i], &sM[i]);
292 // matlab used indexing from 1, C++ starts at 0
293 // so subtract 1:
294 iM[i]--;
295 jM[i]--;
296 // cout << "iM: " << iM[i]
297 // << "\tjM: " << jM[i]
298 // << "\tsM: " << sM[i]
299 // << endl;
300 }
301
303 int offset = 0;
304 for (i=0; i<NumGlobalElements; i++)
305 {
306 // cout << "NumNz[" << i << "] = " << NumNz[i] << endl;
307 for (j=0; j<NumNz[i]; j++)
308 {
310 Indices[j] = jM[offset + j];
311 Values[j] = sM[offset + j];
312 // cout << "iM: " << iM[offset + j]
313 // << "\tjM: " << jM[offset + j]
314 // << "\tsM: " << sM[offset + j]
315 // << endl;
316 }
317 NumEntries = NumNz[i];
318 assert(M.InsertGlobalValues(MyGlobalElements[i],
319 NumEntries, Values, Indices)==0);
320 // cout << "offset = " << offset << endl;
321 offset = offset + NumNz[i];
322 // cout << "Got to here in B matrix." <<endl;
323 }
324
325 // Finish up
326 assert(M.FillComplete(ColMap, Map)==0);
327
328 cout << "nonzeros = " << M.NumGlobalNonzeros() << endl;
329 cout << "rows = " << M.NumGlobalRows() << endl;
330 cout << "cols = " << M.NumGlobalCols() << endl;
331 return Mptr;
332}
333
334
335
336
337
338Petra_RDP_Vector* readVectorIn(FILE *dataFile, Petra_Comm& Comm)
339{
345
346 int NumGlobalElements = 0;
347 int NzElms = 0;
348 int i;
349
350 fscanf(dataFile, "%d", &NumGlobalElements);
351 fscanf(dataFile, "%d", &NzElms);
352 // Construct a Map that puts approximately the same number of
353 // equations on each processor.
354 Petra_Map& Map = *new Petra_Map(NumGlobalElements, 0, Comm);
355
356 // Get update list and number of local equations from newly created Map.
357 int NumMyElements = Map.NumMyElements();
358 int * MyGlobalElements = new int[NumMyElements];
359 Map.MyGlobalElements(MyGlobalElements);
360
361 // make a petra map filled with zeros
362 Petra_RDP_Vector* vptr = new Petra_RDP_Vector(Map);
363 Petra_RDP_Vector& v = *vptr;
364
365 // now fill in the nonzero elements
366 double * myArray = new double[NumMyElements];
367 int tempInd;
368 double tempVal;
369 // cout << "Length v " << NumGlobalElements << endl;
370 // cout << "NzElms " << NzElms << endl;
371 for (i=0; i<NzElms; i++)
372 {
373 fscanf(dataFile, "%d %lg", &tempInd, &tempVal);
374 v[tempInd] = tempVal;
375 // cout << tempVal << endl;
376 }
377 // Petra_RDP_CRS_Matrix& M = *new Petra_RDP_CRS_Matrix(Copy, Map, NumNz);
378
379 return vptr;
380
381}
382
383// small matrix vector multiply test
384void matVecTest(Petra_RDP_CRS_Matrix* Aptr,
385 Petra_RDP_Vector* xptr,
386 Petra_RDP_Vector* bptr)
387{
388 Petra_RDP_CRS_Matrix& A = *Aptr;
389 Petra_RDP_Vector& x = *xptr;
390 Petra_RDP_Vector& b = *bptr; // we're going to overwrite b anyway
391 // will that overwrite x, too? Look at ExtractCopy for alternative.
392
393 A.Multiply(1, x, b);
394}
395
396
397
398// matrix vector multiply for our block matrix
399/************
400Petra_RDP_Vector* = myMatVecMult(Petra_RDP_CRS_Matrix* Bptr,
401 Petra_RDP_CRS_Matrix* Fptr,
402 Petra_RDP_CRS_Matrix* Cptr,
403 Petra_RDP_Vector* xptr)
404
405{
406 // A = [F B'; B C]
407 // return Ax pointer
408 // cout << "block matrix-vector multiply" << endl;
409
410}
411*************/
@ Copy
Petra_RDP_CRS_Matrix * readMatrixIn(FILE *dataFile, Petra_Comm &Comm)
int main(int argc, char *argv[])
Petra_RDP_CRS_Matrix * readRectMatrixIn(FILE *dataFile, Petra_Comm &Comm)
void matVecTest(Petra_RDP_CRS_Matrix *Aptr, Petra_RDP_Vector *xptr, Petra_RDP_Vector *bptr)
Petra_RDP_Vector * readVectorIn(FILE *dataFile, Petra_Comm &Comm)