Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
cc_main.cc
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25//
26// ***********************************************************************
27//@HEADER
28*/
29
30#include <stdio.h>
31#include <stdlib.h>
32#include <ctype.h>
33#include <assert.h>
34#include <string.h>
35#include <math.h>
36#ifdef PETRA_MPI
37#include "mpi.h"
38#endif
39#include "prototypes.h"
40#include "paz_aztec.h"
41#ifndef __cplusplus
42#define __cplusplus
43#endif
44#include "Petra_Comm.h"
45#include "Petra_Map.h"
46#include "Petra_Time.h"
47#include "Petra_BlockMap.h"
48#include "Petra_RDP_MultiVector.h"
49#include "Petra_RDP_Vector.h"
50#include "Petra_RDP_DVBR_Matrix.h"
51#include "Petra_RDP_CRS_Matrix.h"
52#include "Ifpack_ILUK_Graph.h"
53#include "Ifpack_RDP_CRS_RILUK.h"
54
55#define perror(str) { fprintf(stderr,"%s\n",str); exit(-1); }
56#define perror1(str,ierr) { fprintf(stderr,"%s %d\n",str,ierr); exit(-1); }
57#define double_quote '"'
58void BiCGSTAB(Petra_RDP_CRS_Matrix &A,
59 Petra_RDP_Vector &x,
60 Petra_RDP_Vector &b,
61 Ifpack_RDP_CRS_RILUK *M,
62 int Maxiter,
63 double Tolerance,
64 double *residual, double & FLOPS, bool verbose);
65
66int main(int argc, char *argv[])
67{
68 using std::cout;
69 using std::endl;
70
71 int *update; /* vector elements updated on this node. */
72 int *indx; /* MSR format of real and imag parts */
73 int *bindx;
74 int *bpntr;
75 int *rpntr;
76 int *cpntr;
77 int indexBase = 0;
78 double *val;
79 double *xguess, *b, *bt, *xexact, *xsolve;
80 int n_nonzeros, n_blk_nonzeros, ierr;
81 int N_update; /* # of block unknowns updated on this node */
82 int numLocalEquations;
83 /* Number scalar equations on this node */
84 int numGlobalEquations, numGlobalBlocks; /* Total number of equations */
85 int numLocalBlocks;
86 int *blockSizes, *numNzBlks, *blkColInds;
87 int *numNz, *ColInds;
88 int N_external, N_blk_eqns;
89 int blk_row, *blk_col_inds;
90 int row, *col_inds, numEntries;
91 double *row_vals;
92
93 double *val_msr;
94 int *bindx_msr;
95
96 double norm, d ;
97
98 int matrix_type;
99
100 int has_global_indices, option;
101 int i, j, m, mp;
102 int ione = 1;
103 int *proc_config = new int[PAZ_PROC_SIZE];
104
105 double time ;
106#ifdef PETRA_MPI
107 MPI_Init(&argc,&argv);
108#endif
109
110 /* get number of processors and the name of this processor */
111
112#ifdef PETRA_MPI
113 Petra_Comm& Comm = *new Petra_Comm(MPI_COMM_WORLD);
114#else
115 Petra_Comm& Comm = *new Petra_Comm();
116#endif
117
118 printf("proc %d of %d is alive\n",
119 Comm.MyPID(),Comm.NumProc());
120
121 int MyPID = Comm.MyPID();
122
123 bool verbose = false;
124 if (MyPID==0) verbose = true;
125
126
127/* Still need Aztec proc info for the HB routines, can switch to Petra later */
128
129#ifdef PETRA_MPI
130 PAZ_set_proc_config(proc_config,MPI_COMM_WORLD);
131#else
132 PAZ_set_proc_config(proc_config,0);
133#endif
134
135 Comm.Barrier();
136
137#ifdef VBRMATRIX
138 if(argc != 6)
139 perror("error: enter name of data and partition file on command line, followed by levelfill and shift value") ;
140#else
141 if(argc != 5) perror("error: enter name of data file on command line, followed by levelfill and shift value") ;
142#endif
143 /* Set exact solution to NULL */
144 xexact = NULL;
145
146 /* Read matrix file and distribute among processors.
147 Returns with this processor's set of rows */
148
149#ifdef VBRMATRIX
150 if (Comm.MyPID()==0)
151 read_hb(argv[1], &numGlobalEquations, &n_nonzeros,
152 &val_msr, &bindx_msr, &xguess, &b, &bt, &xexact);
153
154 create_vbr(argv[2], proc_config, &numGlobalEquations, &numGlobalBlocks,
155 &n_nonzeros, &n_blk_nonzeros, &N_update, &update,
156 bindx_msr, val_msr, &val, &indx,
157 &rpntr, &cpntr, &bpntr, &bindx);
158
159 if(proc_config[PAZ_node] == 0)
160 {
161 free ((void *) val_msr);
162 free ((void *) bindx_msr);
163 free ((void *) cpntr);
164 }
165 matrix_type = PAZ_VBR_MATRIX;
166
167 Comm.Barrier();
168
169 distrib_vbr_matrix( proc_config, numGlobalEquations, numGlobalBlocks,
170 &n_nonzeros, &n_blk_nonzeros,
171 &N_update, &update,
172 &val, &indx, &rpntr, &cpntr, &bpntr, &bindx,
173 &xguess, &b, &bt, &xexact);
174
175#else
176 if (Comm.MyPID()==0)
177 read_hb(argv[1], &numGlobalEquations, &n_nonzeros,
178 &val, &bindx, &xguess, &b, &bt, &xexact);
179
180 Comm.Barrier();
181
182 distrib_msr_matrix(proc_config, &numGlobalEquations, &n_nonzeros, &N_update,
183 &update, &val, &bindx, &xguess, &b, &bt, &xexact);
184
185#ifdef DEBUG
186 for (i = 0; i<N_update; i++)
187 if (val[i] == 0.0 ) printf("Zero diagonal at row %d\n",i);
188#endif
189 matrix_type = PAZ_MSR_MATRIX;
190#endif
191
192#ifdef VBRMATRIX
193 numLocalEquations = rpntr[N_update];
194#else
195 numLocalEquations = N_update;
196#endif
197
198#ifdef VBRMATRIX
199
200 /******** Make integer data structures needed for this interface *********/
201
202 /* Make blockSizes - number of equations in each block row on this proc */
203
204 numLocalBlocks = N_update;
205
206 blockSizes = new int[numLocalBlocks];
207 for (i=0; i<numLocalBlocks; i++)
208 blockSizes[i] = rpntr[i+1]-rpntr[i];
209
210 /* Make numNzBlks - number of block entries in each block row */
211
212 numNzBlks = new int[numLocalBlocks];
213 for (i=0; i<numLocalBlocks; i++)
214 numNzBlks[i] = bpntr[i+1] - bpntr[i];
215
216 /* Make blkColInds - Exactly bindx (just copy pointer) */
217 blkColInds = bindx;
218
219
220 Petra_BlockMap& map = *new Petra_BlockMap(numGlobalEquations, numLocalEquations,
221 update, indexBase, Comm, numGlobalBlocks, numLocalBlocks,
222 blockSizes);
223
224 Petra_RDP_DVBR_Matrix& A = *new Petra_RDP_DVBR_Matrix(map);
225
226 if ((ierr=A.allocate(numNzBlks, blkColInds)))
227 perror1("Error in DVBR_Matrix_allocate:",ierr);
228
229 /* Add block rows one-at-a-time */
230
231 for (blk_row=0; blk_row<numLocalBlocks; blk_row++)
232 {
233 row_vals = val + indx[bpntr[blk_row]];
234 blk_col_inds = bindx + bpntr[blk_row];
235 if ((ierr=A.putBlockRow(update[blk_row], numNzBlks[blk_row], row_vals,
236 blk_col_inds)))
237 {
238 printf("Row %d:",update[row]);
239 perror1("Error putting block row:",ierr);
240 }
241 }
242
243 if ((ierr=A.FillComplete()))
244 perror1("Error in DVBR_Matrix_FillComplete:",ierr);
245
246#else
247 /* Make numNzBlks - number of block entries in each block row */
248
249 numNz = new int[numLocalEquations];
250 for (i=0; i<numLocalEquations; i++) numNz[i] = bindx[i+1] - bindx[i];
251
252 /* Make ColInds - Exactly bindx, offset by diag (just copy pointer) */
253 ColInds = bindx+numLocalEquations+1;
254
255 Petra_Map& map = *new Petra_Map(numGlobalEquations, numLocalEquations,
256 update, 0, Comm);
257
258 Comm.Barrier();
259 Petra_Time & FillTimer = *new Petra_Time(Comm);
260 Petra_RDP_CRS_Matrix& A = *new Petra_RDP_CRS_Matrix(Copy, map, numNz);
261
262 /* Add rows one-at-a-time */
263
264 for (row=0; row<numLocalEquations; row++)
265 {
266 row_vals = val + bindx[row];
267 col_inds = bindx + bindx[row];
268 numEntries = bindx[row+1] - bindx[row];
269 if ((ierr = A.InsertGlobalValues(update[row], numEntries, row_vals, col_inds)))
270 {
271 printf("Row %d:",update[row]);
272 perror1("Error putting row:",ierr);
273 }
274 if ((ierr=(A.InsertGlobalValues(update[row], 1, val+row, update+row)<0)))
275 perror1("Error putting diagonal",ierr);
276 }
277 Comm.Barrier();
278 double FillTime = FillTimer.ElapsedTime();
279 if ((ierr=A.FillComplete()))
280 perror1("Error in Petra_RDP_Vector_fillComplete",ierr);
281 double FillCompleteTime = FillTimer.ElapsedTime() - FillTime;
282 if (Comm.MyPID()==0) {
283 cout << "\n\n****************************************************" << endl;
284 cout << "\n\nMatrix construction time (sec) = " << FillTime<< endl;
285 cout << " Matrix FillComplete time (sec) = " << FillCompleteTime << endl;
286 cout << " Total construction time (sec) = " << FillTime+FillCompleteTime << endl<< endl;
287 }
288
289
290
291#endif
292
293#ifdef MULTI_VECTOR
294
295 // Make a second x and b vector that are 2x original x and b; make into a 2-vector.
296
297 int nrhs = 2;
298 double **allx = new double*[nrhs];
299 double *xexact2 = new double[numLocalEquations];
300 for (i = 0;i<numLocalEquations; i++) xexact2[i] = 2.0 * xexact[i];
301 allx[0] = xexact; allx[1] = xexact2;
302
303 double **allb = new double*[nrhs];
304 double *b2 = new double[numLocalEquations];
305 for (i = 0;i<numLocalEquations; i++) b2[i] = 2.0 * b[i];
306 allb[0] = b; allb[1] = b2;
307
308 double **allbt = new double*[nrhs];
309 double *bt2 = new double[numLocalEquations];
310 for (i = 0;i<numLocalEquations; i++) bt2[i] = 2.0 * bt[i];
311 allbt[0] = bt; allbt[1] = bt2;
312
313 Petra_RDP_MultiVector& xtrue = *new Petra_RDP_MultiVector(Copy, map, allx, nrhs);
314 Petra_RDP_MultiVector& bexact = *new Petra_RDP_MultiVector(Copy, map, allb, nrhs);
315 Petra_RDP_MultiVector& btexact = *new Petra_RDP_MultiVector(Copy, map, allbt, nrhs);
316 Petra_RDP_MultiVector& bcomp = *new Petra_RDP_MultiVector(map, nrhs);
317 Petra_RDP_MultiVector& xcomp = *new Petra_RDP_MultiVector(map, nrhs);
318 Petra_RDP_MultiVector& resid = *new Petra_RDP_MultiVector(map, nrhs);
319
320#else
321 int nrhs = 1;
322 Petra_RDP_Vector& xtrue = *new Petra_RDP_Vector(Copy, map, xexact);
323 Petra_RDP_Vector& bexact = *new Petra_RDP_Vector(Copy, map, b);
324 Petra_RDP_Vector& btexact = *new Petra_RDP_Vector(Copy, map, bt);
325 Petra_RDP_Vector& bcomp = *new Petra_RDP_Vector(map);
326 Petra_RDP_Vector& xcomp = *new Petra_RDP_Vector(map);
327 Petra_RDP_Vector& resid = *new Petra_RDP_Vector(map);
328
329#endif /* MULTI_VECTOR */
330
331 Comm.Barrier();
332
333 // Construct ILU preconditioner
334
335 double elapsed_time, total_flops, MFLOPs;
336 Petra_Time & timer = *new Petra_Time(Comm);
337
338 int LevelFill = atoi(argv[argc-3]);
339 if (verbose) cout << "Using Level Fill = " << LevelFill << endl;
340 double ShiftValue = atof(argv[argc-2]);
341 if (verbose) cout << "Using Diagonal Shift Value of = " << ShiftValue << endl;
342 int NumThreads = atoi(argv[argc-1]);
343 if (verbose) cout << "Using " << NumThreads << " Threads " << endl;
344
345 Ifpack_RDP_CRS_RILUK * ILUK = 0;
346 Ifpack_ILUK_Graph * ILUK_Graph = 0;
347 if (LevelFill>-1) {
348 elapsed_time = timer.ElapsedTime();
349 ILUK_Graph = new Ifpack_ILUK_Graph(A.Graph(), LevelFill);
350 assert(ILUK_Graph->ConstructFilledGraph()==0);
351 assert(ILUK_Graph->ComputeLevels(NumThreads)==0);
352 elapsed_time = timer.ElapsedTime() - elapsed_time;
353 if (verbose) cout << "Time to construct ILUK graph = " << elapsed_time << endl;
354
355 return 0;
356
357 elapsed_time = timer.ElapsedTime();
358 ILUK = new Ifpack_RDP_CRS_RILUK(A, *ILUK_Graph);
359 ILUK->SetShiftValue(ShiftValue);
360 assert(ILUK->InitValues()==0);
361 assert(ILUK->Factor()==0);
362 elapsed_time = timer.ElapsedTime() - elapsed_time;
363 total_flops = ILUK->Flops();
364 MFLOPs = total_flops/elapsed_time/1000000.0;
365 if (verbose) cout << "Time to compute preconditioner values = "
366 << elapsed_time << endl
367 << "MFLOPS for Factorization = " << MFLOPs << endl;
368 }
369
370 int Maxiter = 500;
371 double Tolerance = 1.0E-14;
372 double * residual = new double[nrhs];
373
374
375 elapsed_time = timer.ElapsedTime();
376
377 BiCGSTAB(A, xcomp, bexact, ILUK, Maxiter, Tolerance, residual, total_flops, verbose);
378
379 elapsed_time = timer.ElapsedTime() - elapsed_time;
380 MFLOPs = total_flops/elapsed_time/1000000.0;
381 if (verbose) cout << "Time to compute solution = "
382 << elapsed_time << endl
383 << "Number of operations in solve = " << total_flops << endl
384 << "MFLOPS for Solve = " << MFLOPs<< endl << endl;
385
386
387 int NumMVs = 100;
388 int ii, iii;
389
390 for (ii=0; ii<2; ii++) {
391 bool TransA = (ii==1);
392 A.ResetFlops();
393 elapsed_time = timer.ElapsedTime();
394 for (iii=0; iii<NumMVs; iii++)
395 if ((ierr=A.Multiply(TransA, xcomp, bcomp))) perror1("Error in matvec",ierr);
396 elapsed_time = timer.ElapsedTime() - elapsed_time;
397 total_flops = A.Flops();
398 MFLOPs = total_flops/elapsed_time/1000000.0;
399 if (Comm.MyPID()==0) {
400 if (TransA) {
401 cout << "\n\n****************************************************" << endl;
402 cout << "\n\nResults for transpose multiply with standard storage" << endl;
403 }
404 else {
405 cout << "\n\nMatrix Fill cost = " << (FillTime+FillCompleteTime)/elapsed_time*NumMVs
406 << " Matrix-Multiplies " << endl;
407 cout << "\n\n*******************************************************" << endl;
408 cout << "\n\nResults for no transpose multiply with standard storage" << endl;
409 }
410
411 cout << "\n\nMatrix Fill cost = " << (FillTime+FillCompleteTime)/elapsed_time*NumMVs
412 << " Matrix-Multiplies " << endl;
413 cout << "\n\nTotal FLOPS for standard Storage (" <<NumMVs<< " Multiplies) = "
414 << total_flops<< endl;
415 cout << " Total time for standard Storage = " << elapsed_time<< endl;
416 cout << " Total MFLOPs for standard matrix multiply = " << MFLOPs << endl<< endl;
417 }
418
419 // cout << "Vector bcomp" << bcomp << endl;
420
421 if (TransA) {
422 if ((ierr=resid.Update(-1.0, btexact, 1.0, bcomp, 0.0))) perror1("Error in Update",ierr);}
423 else {
424 if ((ierr=resid.Update(-1.0, bexact, 1.0, bcomp, 0.0))) perror1("Error in Update",ierr);}
425
426 if ((ierr = resid.Norm2(residual))) perror1("Error in Norm2",ierr);
427
428 if (Comm.MyPID()==0)
429 for (i=0; i< nrhs; i++) printf("Residual[%d] = %22.16g\n",i,residual[i]);
430
431 }
432
433 // Optimize data layout for memory access
434
435 if ((ierr=A.OptimizeStorage())) perror1("Error in Petra_RDP_CRS_Matrix: OptimizeStorage",ierr);
436
437 for (ii=0; ii<2; ii++) {
438 bool TransA = (ii==1);
439 A.ResetFlops();
440 elapsed_time = timer.ElapsedTime();
441 for (iii=0; iii<NumMVs; iii++)
442 if ((ierr=A.Multiply(TransA, xcomp, bcomp)))
443 perror1("Error in Multiply",ierr);
444 elapsed_time = timer.ElapsedTime() - elapsed_time;
445 total_flops = A.Flops();
446 MFLOPs = total_flops/elapsed_time/1000000.0;
447 if (Comm.MyPID()==0) {
448 cout << "\n\n****************************************************" << endl;
449 if (TransA) cout << "\n\nResults for transpose multiply with optimized storage" << endl;
450 else cout << "\n\nResults for no transpose multiply with optimized storage"<< endl;
451 cout << "\n\nTotal FLOPS for optimized storage (" <<NumMVs<< " Multiplies) = "
452 << total_flops<< endl;
453 cout << " Total time for optimized Storage = " << elapsed_time<< endl;
454 cout << " Total MFLOPs for optimized matrix multiply = " << MFLOPs << endl<< endl;
455 }
456
457 if (TransA) {
458 if ((ierr=resid.Update(-1.0, btexact, 1.0, bcomp, 0.0))) perror1("Error in Update",ierr);}
459 else {
460 if ((ierr=resid.Update(-1.0, bexact, 1.0, bcomp, 0.0))) perror1("Error in Update",ierr); }
461
462 if ((ierr = resid.Norm2(residual))) perror1("Error in Norm2",ierr);
463
464 if (Comm.MyPID()==0)
465 for (i=0; i< nrhs; i++) printf("Residual[%d] = %22.16g\n",i,residual[i]);
466
467 }
468
469 free ((void *) xguess);
470 free ((void *) b);
471 free ((void *) bt);
472 free ((void *) xexact);
473 free ((void *) val);
474 free ((void *) bindx);
475 free ((void *) update);
476
477#ifdef VBRMATRIX
478 free ((void *) indx);
479 free ((void *) rpntr);
480 free ((void *) bpntr);
481 delete [] blockSizes;
482 delete [] numNzBlks;
483#else
484 delete [] numNz;
485
486#endif
487
488#ifdef MULTI_VECTOR
489 delete [] xexact2;
490 delete [] b2;
491 delete [] allx;
492 delete [] allb;
493#endif
494 delete &bexact;
495 delete &bcomp;
496 delete &resid;
497 delete &xcomp;
498 delete ILUK;
499 delete &ILUK_Graph;
500 delete &A;
501 delete &map;
502 delete &timer;
503 delete &FillTimer;
504 delete &Comm;
505
506 delete proc_config;
507
508#ifdef PETRA_MPI
509 MPI_Finalize() ;
510#endif
511
512return 0 ;
513}
514void BiCGSTAB(Petra_RDP_CRS_Matrix &A,
515 Petra_RDP_Vector &x,
516 Petra_RDP_Vector &b,
517 Ifpack_RDP_CRS_RILUK *M,
518 int Maxiter,
519 double Tolerance,
520 double *residual, double & FLOPS, bool verbose) {
521
522 // Allocate vectors needed for iterations
523 Petra_RDP_Vector& phat = *new Petra_RDP_Vector(x); phat.PutScalar(0.0);
524 Petra_RDP_Vector& p = *new Petra_RDP_Vector(x); p.PutScalar(0.0);
525 Petra_RDP_Vector& shat = *new Petra_RDP_Vector(x); shat.PutScalar(0.0);
526 Petra_RDP_Vector& s = *new Petra_RDP_Vector(x); s.PutScalar(0.0);
527 Petra_RDP_Vector& r = *new Petra_RDP_Vector(x); r.PutScalar(0.0);
528 Petra_RDP_Vector& rtilde = *new Petra_RDP_Vector(x); rtilde.Random();
529 Petra_RDP_Vector& v = *new Petra_RDP_Vector(x); v.PutScalar(0.0);
530
531
532 A.Multiply(false, x, r); // r = A*x
533
534 r.Update(1.0, b, -1.0); // r = b - A*x
535
536 double r_norm, b_norm, scaled_r_norm, rhon, rhonm1 = 1.0;
537 double alpha = 1.0, omega = 1.0, sigma;
538 double omega_num, omega_den;
539 r.Norm2(&r_norm);
540 b.Norm2(&b_norm);
541 scaled_r_norm = r_norm/b_norm;
542 r.Dot(rtilde,&rhon);
543
544 if (verbose) cout << "Initial residual = " << r_norm
545 << " Scaled residual = " << scaled_r_norm << endl;
546
547
548 for (int i=0; i<Maxiter; i++) { // Main iteration loop
549
550 double beta = (rhon/rhonm1) * (alpha/omega);
551 rhonm1 = rhon;
552
553 /* p = r + beta*(p - omega*v) */
554 /* phat = M^-1 p */
555 /* v = A phat */
556
557 double dtemp = - beta*omega;
558
559 p.Update(1.0, r, dtemp, v, beta);
560 if (M==0)
561 phat.Scale(1.0, p);
562 else
563 M->LevelSolve(false, p, phat);
564 A.Multiply(false, phat, v);
565
566
567 rtilde.Dot(v,&sigma);
568 alpha = rhon/sigma;
569
570 /* s = r - alpha*v */
571 /* shat = M^-1 s */
572 /* r = A shat (r is a tmp here for t ) */
573
574 s.Update(-alpha, v, 1.0, r, 0.0);
575 if (M==0)
576 shat.Scale(1.0, s);
577 else
578 M->LevelSolve(false, s, shat);
579 A.Multiply(false, shat, r);
580
581 r.Dot(s, &omega_num);
582 r.Dot(r, &omega_den);
583 omega = omega_num / omega_den;
584
585 /* x = x + alpha*phat + omega*shat */
586 /* r = s - omega*r */
587
588 x.Update(alpha, phat, omega, shat, 1.0);
589 r.Update(1.0, s, -omega);
590
591 r.Norm2(&r_norm);
592 scaled_r_norm = r_norm/b_norm;
593 r.Dot(rtilde,&rhon);
594
595 if (verbose) cout << "Iter "<< i << " residual = " << r_norm
596 << " Scaled residual = " << scaled_r_norm << endl;
597
598 if (scaled_r_norm < Tolerance) break;
599 }
600
601 FLOPS = phat.Flops()+p.Flops()+shat.Flops()+s.Flops()+r.Flops()+rtilde.Flops()+
602 v.Flops()+A.Flops()+x.Flops()+b.Flops();
603 if (M!=0) FLOPS += M->Flops();
604
605 delete &phat;
606 delete &p;
607 delete &shat;
608 delete &s;
609 delete &r;
610 delete &rtilde;
611 delete &v;
612
613
614 return;
615}