69int main(
int argc,
char *argv[])
79 MPI_Init(&argc,&argv);
82 MPI_Comm_size(MPI_COMM_WORLD, &size);
83 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
104 if (rank==0) cout <<
"Press any key to continue..."<< endl;
105 if (rank==0) cin >> tmp;
108 int MyPID = Comm.
MyPID();
110 if (verbose) cout << Comm << endl;
112 bool verbose1 = verbose;
115 if (verbose && rank!=0) verbose =
false;
118 if(argc != 2) cout <<
"Error: Enter name of data file on command line." << endl;
123 int NumGlobalEquations, n_nonzeros, *bindx;
124 double * val, * xguess, * b, * xexact = 0;
126 Trilinos_Util_read_hb(argv[1], Comm.
MyPID(), &NumGlobalEquations, &n_nonzeros,
127 &val, &bindx, &xguess, &b, &xexact);
129 int NumMyEquations, * MyGlobalEquations;
131 Trilinos_Util_distrib_msr_matrix(Comm, &NumGlobalEquations, &n_nonzeros, &NumMyEquations,
132 &MyGlobalEquations, &val, &bindx, &xguess, &b, &xexact);
137 int * NumNz =
new int[NumMyEquations];
138 for (i=0; i<NumMyEquations; i++) NumNz[i] = bindx[i+1] - bindx[i] + 1;
141 Epetra_Map Map(NumGlobalEquations, NumMyEquations,
142 MyGlobalEquations, 0, Comm);
154 for (i=0; i<NumMyEquations; i++){
155 Values = val + bindx[i];
156 Indices = bindx + bindx[i];
157 NumIndices = bindx[i+1] - bindx[i];
164 if (verbose) cout <<
"\nTime to construct A = " << timer.
ElapsedTime() - start << endl;
165 double * xexactt = xexact;
173 assert(A.
Multiply(
false, xx, bcomp)==0);
176 assert(resid.
Update(1.0, bb, -1.0, bcomp, 0.0)==0);
179 assert(resid.
Norm2(&residual)==0);
180 if (Comm.
MyPID()==0) cout <<
"Sanity check: Norm of b - Ax for known x and b = " << residual << endl;
195 int * TransNumNz =
new int[NumMyCols];
196 for (i=0;i<NumMyCols; i++) TransNumNz[i] = 0;
197 for (i=0; i<NumMyEquations; i++) {
199 for (j=0; j<NumIndices; j++) ++TransNumNz[Indices[j]];
202 int ** TransIndices =
new int*[NumMyCols];
203 double ** TransValues =
new double*[NumMyCols];
205 for(i=0; i<NumMyCols; i++) {
206 NumIndices = TransNumNz[i];
208 TransIndices[i] =
new int[NumIndices];
209 TransValues[i] =
new double[NumIndices];
215 for (i=0;i<NumMyCols; i++) TransNumNz[i] = 0;
216 for (i=0; i<NumMyEquations; i++) {
219 for (j=0; j<NumIndices; j++) {
220 int TransRow = Indices[j];
221 int loc = TransNumNz[TransRow];
222 TransIndices[TransRow][loc] = ii;
223 TransValues[TransRow][loc] = Values[j];
224 ++TransNumNz[TransRow];
234 int * TransMyGlobalEquations =
new int[NumMyCols];
239 for (i=0; i<NumMyCols; i++)
242 TransNumNz[i], TransValues[i], TransIndices[i])==0);
258 assert(TransA1.
Export(TempTransA1, Export,
Add)==0);
263 if (verbose) cout <<
"\nTime to construct TransA1 = " << timer.
ElapsedTime() - start << endl;
274 assert(TransA1.
Multiply(
false, x, b2)==0);
276 assert(resid.
Update(1.0, b1, -1.0, b2, 0.0)==0);
278 assert(b1.
Norm2(&residual)==0);
279 if (verbose) cout <<
"Norm of RHS using Trans = true with A = " << residual << endl;
280 assert(b2.
Norm2(&residual)==0);
281 if (verbose) cout <<
"Norm of RHS using Trans = false with TransA1 = " << residual << endl;
282 assert(resid.
Norm2(&residual)==0);
283 if (verbose) cout <<
"Difference between using A and TransA1 = " << residual << endl;
299 for (
int LocalRow=0; LocalRow<NumMyEquations; LocalRow++) {
301 int TransGlobalCol = A.
GRID(LocalRow);
302 for (j=0; j<NumIndices; j++) {
303 int TransGlobalRow = A.
GCID(Indices[j]);
304 double TransValue = Values[j];
305 assert(TempTransA2.
InsertGlobalValues(TransGlobalRow, 1, &TransValue, &TransGlobalCol)>=0);
315 if (verbose) cout <<
"\nTime to construct TransA2 = " << timer.
ElapsedTime() - start << endl;
326 assert(TransA2.
Export(TempTransA2, Export,
Add)==0);
339 assert(TransA2.
Multiply(
false, x, b2)==0);
341 assert(resid.
Update(1.0, b1, -1.0, b2, 0.0)==0);
343 assert(b1.
Norm2(&residual)==0);
344 if (verbose) cout <<
"Norm of RHS using Trans = true with A = " << residual << endl;
345 assert(b2.
Norm2(&residual)==0);
346 if (verbose) cout <<
"Norm of RHS using Trans = false with TransA2 = " << residual << endl;
347 assert(resid.
Norm2(&residual)==0);
348 if (verbose) cout <<
"Difference between using A and TransA2 = " << residual << endl;
352 free ((
void *) xguess);
354 free ((
void *) xexact);
356 free ((
void *) bindx);
357 free ((
void *) MyGlobalEquations);
359 delete [] TransMyGlobalEquations;
361 for(i=0; i<NumMyCols; i++) {
362 NumIndices = TransNumNz[i];
364 delete [] TransIndices[i];
365 delete [] TransValues[i];
368 delete [] TransIndices;
369 delete [] TransValues;
372 delete [] TransNumNz;
int main(int argc, char *argv[])