]> AND Private Git Repository - GMRES2stage.git/blob - code/ex10.c
Logo AND Algorithmique Numérique Distribuée

Private GIT Repository
update experiments
[GMRES2stage.git] / code / ex10.c
1
2 ///home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/crashbasis/crashbasis.bin  -ksp_type gmres  -pc_type none 
3 //16s
4 //temps simi pour lui et nous
5
6
7 // /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/parabolic_fem/parabolic_fem.bin  -ksp_type gmres  -pc_type ilu
8 //converge avec nous mais lentement
9 //2152s pour lui
10 //724s pour nous
11
12
13
14 // /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/epb3/epb3.bin  -ksp_type fgmres  -pc_type sor
15 //fonctionne avec ilu asm 
16 //9s
17 //temps simi pour lui et nous
18
19 // /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/atmosmodj/atmosmodj.bin  -ksp_type fgmres  -pc_type none
20 //aussi sor
21
22 // /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/bfwa398/bfwa398.bin  -ksp_type gmres   -pc_type none
23
24
25 // /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 1 ./ex10 -f ~/Downloads/torso3/torso3.bin  -ksp_type gmres  -pc_type none
26
27
28 static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
29 This version first preloads and solves a small system, then loads \n\
30 another (larger) system and solves it as well.  This example illustrates\n\
31 preloading of instructions with the smaller system so that more accurate\n\
32 performance monitoring can be done with the larger one (that actually\n\
33 is the system of interest).  See the 'Performance Hints' chapter of the\n\
34 users manual for a discussion of preloading.  Input parameters include\n\
35   -f0 <input_file> : first file to load (small system)\n\
36   -f1 <input_file> : second file to load (larger system)\n\n\
37   -nearnulldim <0> : number of vectors in the near-null space immediately following matrix\n\n\
38   -trans  : solve transpose system instead\n\n";
39 /*
40   This code can be used to test PETSc interface to other packages.\n\
41   Examples of command line options:       \n\
42    ./ex10 -f0 <datafile> -ksp_type preonly  \n\
43         -help -ksp_view                  \n\
44         -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
45         -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu or superlu_dist or mumps \n\
46         -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package mumps \n\
47    mpiexec -n <np> ./ex10 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
48  \n\n";
49 */
50 /*T
51    Concepts: KSP^solving a linear system
52    Processors: n
53 T*/
54
55 /*
56   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
57   automatically includes:
58      petscsys.h       - base PETSc routines   petscvec.h - vectors
59      petscmat.h - matrices
60      petscis.h     - index sets            petscksp.h - Krylov subspace methods
61      petscviewer.h - viewers               petscpc.h  - preconditioners
62 */
63 #include <petscksp.h>
64
65 #undef __FUNCT__
66 #define __FUNCT__ "main"
67
68
69
70
71
72
73
74
75
76 int KrylovMinimize(Mat A, Vec b, Vec x) {
77   
78
79   //Variables
80
81   PetscScalar  gamma, alpha, oldgamma, beta;
82   PetscReal norm=20, Eprecision=1e-7, cgprec=1e-40;     
83   PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;
84   PetscErrorCode ierr;
85   PetscScalar T1, T2;
86   KSP ksp;
87   PetscInt total=0;  
88   PetscInt size;
89   PetscInt Istart,Iend;
90   PetscInt i,its;
91   Vec       x_old, residu;
92   Mat S, AS;
93   PetscScalar *array;
94   PetscInt *ind_row;
95   Vec Alpha, p, ss, vect, r, q, Ax;
96
97
98   PetscInt first=1;
99
100   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
101   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
102
103
104
105   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
106
107
108
109
110   VecGetSize(b,&size);
111
112   ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);  
113
114   PetscInt aa,bb;
115   MatGetOwnershipRange(A,&aa,&bb);
116
117   //  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);  
118   //PetscSynchronizedFlush(PETSC_COMM_WORLD);
119
120
121   ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
122   ierr = MatSetSizes(S, bb-aa,  PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
123   ierr = MatSetType(S, MATDENSE); CHKERRQ(ierr);
124   ierr = MatSetUp(S); CHKERRQ(ierr);
125
126   ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
127
128   ierr = VecCreate(PETSC_COMM_WORLD, &Alpha); CHKERRQ(ierr);
129   ierr = VecSetSizes(Alpha, PETSC_DECIDE, ColS); CHKERRQ(ierr);
130   ierr = VecSetFromOptions(Alpha); CHKERRQ(ierr);
131   ierr = VecDuplicate(Alpha, &vect); CHKERRQ(ierr);
132   ierr = VecDuplicate(Alpha, &p); CHKERRQ(ierr);
133   ierr = VecDuplicate(Alpha, &ss); CHKERRQ(ierr);
134   ierr = VecDuplicate(b, &r); CHKERRQ(ierr);
135   ierr = VecDuplicate(b, &q); CHKERRQ(ierr);
136   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
137
138   ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
139   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
140
141
142
143   //indexes of row (these indexes are global)
144   ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
145   for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
146
147   //Initializations
148   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
149   ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, PETSC_DEFAULT, 30); CHKERRQ(ierr);
150   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
151
152
153
154   //GMRES WITH MINIMIZATION
155   T1 = MPI_Wtime();
156   ierr = KSPSetUp(ksp);CHKERRQ(ierr);
157   while(giter<Emaxiter && norm>Eprecision ){
158     for(col=0; col<ColS  &&  norm>Eprecision; col++){
159
160
161       //Solve
162       ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
163
164       ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
165       total += its; Iiter ++;
166
167
168
169       //Build S'
170       ierr = VecGetArray(x, &array);
171       ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
172       VecRestoreArray(x, &array);
173
174
175
176       KSPGetResidualNorm(ksp,&norm);
177
178       /*      //Error
179       ierr = VecCopy(x, residu); CHKERRQ(ierr);
180       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
181       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
182        */
183
184
185       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
186       //ierr = VecCopy(x, x_old); CHKERRQ(ierr);
187
188
189     }
190
191
192     //minimization step
193     if( norm>Eprecision) {
194
195       MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
196       MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
197
198
199       //Build AS
200       if(first) {
201         MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
202         first=0;
203       }
204       else
205         MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
206
207
208
209
210       //Minimization with CGLS  
211       MatMult(AS, Alpha, r);  VecAYPX(r, -1, b); //r_0 = b-AS*x_0
212
213
214       MatMultTranspose(AS, r, p); //p_0 = AS'*r_0 
215
216
217
218
219       ierr = VecCopy(p, ss); CHKERRQ(ierr); //p_0 = s_0
220       VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
221       while(gamma>cgprec && iter<iterations){
222         MatMult(AS, p, q); //q = AS*p
223         VecNorm(q, NORM_2, &alpha); alpha = alpha * alpha; //alpha = norm2(q)^2
224         alpha = gamma / alpha; 
225         VecAXPY(Alpha, alpha, p); //Alpha += alpha*p
226         VecAXPY(r, -alpha, q); //r -= alpha*q
227         MatMultTranspose(AS, r, ss); // ss = AS'*r
228         oldgamma = gamma;
229         VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
230         beta = gamma / oldgamma;
231         VecAYPX(p, beta, ss); //p = s+beta*p;
232         iter ++;
233       } 
234
235       iter = 0; giter ++;
236       //Minimizer
237       MatMult(S, Alpha, x); //x = S*Alpha;
238     }
239
240   }
241   T2 = MPI_Wtime();
242   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
243   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
244
245   return 0;
246
247 }
248
249
250
251
252
253
254
255
256
257 int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
258   
259
260   //Variables
261
262   PetscScalar  alpha, beta;
263   PetscReal norm=20, Eprecision=1e-7, tol=1e-40;     
264   PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=20, Iiter=0;
265   PetscErrorCode ierr;
266   PetscScalar T1, T2;
267   KSP ksp;
268   PetscInt total=0;  
269   PetscInt size;
270   PetscInt Istart,Iend;
271   PetscInt i,its;
272   Vec       x_old, residu;
273   Mat S, AS;
274   PetscScalar *array;
275   PetscInt *ind_row;
276   Vec  Ax;
277   PetscScalar norm_ksp;
278   Vec u,v,d,uu,vt,zero_long,zero_short,x_lsqr;
279
280   PetscInt first=1;
281
282   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
283   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
284   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
285
286
287
288
289   VecGetSize(b,&size);
290
291   ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);  
292
293   PetscInt aa,bb;
294   MatGetOwnershipRange(A,&aa,&bb);
295
296   //  ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%D %D\n", aa,bb); CHKERRQ(ierr);  
297   //PetscSynchronizedFlush(PETSC_COMM_WORLD);
298
299
300   ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
301   ierr = MatSetSizes(S, bb-aa,  PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
302   ierr = MatSetType(S, MATDENSE); CHKERRQ(ierr);
303   ierr = MatSetUp(S); CHKERRQ(ierr);
304
305   ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
306
307
308
309   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
310
311   ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
312   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
313
314
315   //long vector
316   ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
317
318
319   ierr = VecDuplicate(b,&uu);CHKERRQ(ierr);
320   ierr = VecDuplicate(b,&zero_long);CHKERRQ(ierr);
321   ierr = VecSet(zero_long,0);CHKERRQ(ierr);
322
323   //small vector
324   ierr = VecCreate(PETSC_COMM_WORLD, &v); CHKERRQ(ierr);
325   ierr = VecSetSizes(v, PETSC_DECIDE, ColS); CHKERRQ(ierr);
326   ierr = VecSetFromOptions(v); CHKERRQ(ierr);
327   ierr = VecDuplicate(v,&zero_short);CHKERRQ(ierr);
328   ierr = VecSet(zero_short,0);CHKERRQ(ierr);
329   ierr = VecDuplicate(v,&d);CHKERRQ(ierr);
330   ierr = VecDuplicate(v,&vt);CHKERRQ(ierr);
331   ierr = VecDuplicate(v,&x_lsqr);CHKERRQ(ierr);
332
333
334   //indexes of row (these indexes are global)
335   ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
336   for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
337
338   //Initializations
339   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
340   ierr = KSPSetTolerances(ksp, 1e-10, 1e-10, PETSC_DEFAULT, 30); CHKERRQ(ierr);
341   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
342
343
344
345
346   //GMRES WITH MINIMIZATION
347   T1 = MPI_Wtime();
348   ierr = KSPSetUp(ksp);CHKERRQ(ierr);
349   while(giter<Emaxiter && norm>Eprecision ){
350     for(col=0; col<ColS  &&  norm>Eprecision; col++){
351
352
353       //Solve
354       ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
355
356       ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
357       total += its; Iiter ++;
358
359
360
361       //Build S'
362       ierr = VecGetArray(x, &array);
363       ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
364       VecRestoreArray(x, &array);
365
366
367
368       KSPGetResidualNorm(ksp,&norm);
369
370       /*
371       //Error
372       ierr = VecCopy(x, residu); CHKERRQ(ierr);
373       ierr = VecAXPY(residu, -1, x_old); CHKERRQ(ierr);
374       ierr = VecNorm(residu, NORM_INFINITY, &norm); CHKERRQ(ierr);
375        */
376
377
378       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
379       //ierr = VecCopy(x, x_old); CHKERRQ(ierr);
380
381
382     }
383
384
385     //minimization step
386     if( norm>Eprecision) {
387
388       MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
389       MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
390
391
392
393
394       //Build AS
395       if(first) {
396         MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
397         
398         first=0;
399       }
400       else
401         MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
402
403
404
405
406
407       //LSQR
408       //LSQR
409       //LSQR
410
411
412
413       PetscScalar n2b,tolb,normr,c,s,phibar,normar,norma,thet,rhot,rho,phi;
414       PetscInt stag;
415       tolb = tol * n2b;
416       VecNorm(b, NORM_2, &n2b); //n2b = norm(b);
417       ierr = VecCopy(b, u); CHKERRQ(ierr);   //u=b
418       VecNorm(u, NORM_2, &beta); // beta=norm(u)
419       normr=beta;
420       if (beta != 0) {
421         VecAYPX(u,1/beta,zero_long); //  u = u / beta;
422       }
423       c=1;
424       s=0;
425       phibar=beta;
426       MatMultTranspose(AS, u, v); //v=A'*u
427       ierr = VecSet(x_lsqr,0);CHKERRQ(ierr);
428       VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
429       if (alpha != 0) {
430         VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
431       }
432       ierr = VecSet(d,0);CHKERRQ(ierr);
433       normar = alpha * beta;
434       norma=0;
435       //stag=0;
436       for(i=0;i<iterations;i++) {
437         MatMult(AS, v, uu); //uu=A*v
438         VecAYPX(u, -alpha, uu); //u = uu-alpha*u;
439         VecNorm(u, NORM_2, &beta); // beta=norm(u)
440         VecAYPX(u,1/beta,zero_long); //  u = u / beta;
441         norma=sqrt(norma*norma+alpha*alpha+beta*beta); // norma = norm([norma alpha beta]);
442         thet = - s * alpha;
443         rhot = c * alpha;
444         rho = sqrt(rhot*rhot + beta*beta);
445         c = rhot / rho;
446         s = - beta / rho;
447         phi = c * phibar;
448         if (phi == 0) {             // stagnation of the method
449           stag = 1;
450         }
451         phibar = s * phibar;
452         VecAYPX(d,-thet,v);       //d = (v - thet * d);
453         VecAYPX(d,1/rho,zero_short);     //d=d/ rho;
454
455         /*
456         if (normar/(norma*normr) <= tol) { // check for convergence in min{|b-A*x|}
457           break;
458         }
459         if (normr <= tolb) {           // check for convergence in A*x=b
460           break;
461         }
462          */
463
464         VecAXPY(x_lsqr,phi,d);     // x_lsqr=x_lsqr+phi*d
465         normr = abs(s) * normr;
466         MatMultTranspose(AS, u, vt); //vt=A'*u;
467         VecAYPX(v,-beta,vt);  // v = vt - beta * v;
468         VecNorm(v, NORM_2, &alpha); // alpha=norm(v)
469         VecAYPX(v,1/alpha,zero_short); //  v = v / alpha;
470         normar = alpha * abs( s * phi);
471       } 
472       
473
474
475       iter = 0; giter ++;
476       //Minimizer
477       MatMult(S, x_lsqr, x); //x = S*x_small;
478     }
479
480   }
481   T2 = MPI_Wtime();
482   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time LSQR : %g (s)\n", T2-T1); CHKERRQ(ierr);
483   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations LSQR : %D\n", total); CHKERRQ(ierr);
484
485   return 0;
486
487 }
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502 int main(int argc,char **args)
503 {
504   KSP            ksp;             /* linear solver context */
505   Mat            A,A2;               /* matrix */
506   Vec            x,b,u;           /* approx solution, RHS, exact solution */
507   PetscViewer    fd;              /* viewer */
508   char           file[4][PETSC_MAX_PATH_LEN];     /* input file name */
509   PetscBool      table     =PETSC_FALSE,flg,trans=PETSC_FALSE,initialguess = PETSC_FALSE;
510   PetscBool      outputSoln=PETSC_FALSE;
511   PetscErrorCode ierr;
512   PetscInt       its,num_numfac,m,n,M,nearnulldim = 0;
513   PetscReal      norm;
514   PetscBool      preload=PETSC_TRUE,isSymmetric,cknorm=PETSC_FALSE,initialguessfile = PETSC_FALSE;
515   PetscMPIInt    rank;
516   char           initialguessfilename[PETSC_MAX_PATH_LEN];
517
518   PetscInitialize(&argc,&args,(char*)0,help);
519   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
520
521   ierr = PetscOptionsGetString(NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
522   if (flg) {
523     ierr    = PetscStrcpy(file[1],file[0]);CHKERRQ(ierr);
524     preload = PETSC_FALSE;
525   } else {
526     ierr = PetscOptionsGetString(NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
527     if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
528     ierr = PetscOptionsGetString(NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
529     if (!flg) preload = PETSC_FALSE;   /* don't bother with second system */
530   }
531
532
533   PetscPreLoadBegin(preload,"Load system");
534
535
536
537   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);CHKERRQ(ierr);
538
539   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
540   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
541   ierr = MatLoad(A,fd);CHKERRQ(ierr);
542
543   /*  char           ordering[256] = MATORDERINGRCM;
544   IS             rowperm       = NULL,colperm = NULL;
545   ierr = MatGetOrdering(A2,ordering,&rowperm,&colperm);CHKERRQ(ierr);
546   ierr = MatPermute(A2,rowperm,colperm,&A);CHKERRQ(ierr);
547    */
548
549   ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
550
551   ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
552
553  ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
554       ierr = VecSetSizes(b,m,PETSC_DECIDE);CHKERRQ(ierr);
555       ierr = VecSetFromOptions(b);CHKERRQ(ierr);
556       ierr = VecSet(b,1);CHKERRQ(ierr);
557
558
559
560
561   ierr = MatGetVecs(A,&x,NULL);CHKERRQ(ierr);
562   ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
563
564     ierr = VecSet(x,1.0);CHKERRQ(ierr);
565
566
567   ierr       = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
568   ierr       = KSPSetInitialGuessNonzero(ksp,initialguess);CHKERRQ(ierr);
569   num_numfac = 1;
570
571
572
573
574
575     ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
576
577    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
578
579       ierr = KSPSetTolerances(ksp,1e-10,1e-10,PETSC_DEFAULT,5000000);CHKERRQ(ierr);
580   PC             pc;
581       KSPGetPC(ksp, &pc);
582   PCType         type;
583         PCGetType(pc, &type);
584
585          PetscPrintf(PETSC_COMM_WORLD, "PC TYPE %s  \n", type);
586
587
588        PetscScalar T1,T2;
589         T1 = MPI_Wtime();
590         ierr = KSPSetUp(ksp);CHKERRQ(ierr);
591         ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
592         T2 = MPI_Wtime();
593
594         
595     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
596       ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
597
598
599   
600       ierr = MatMult(A,x,u);CHKERRQ(ierr);
601    ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr);
602     ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
603  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr);
604
605  {
606
607     Vec x2;
608     Vec sol;
609     VecDuplicate(b,&x2);
610     VecDuplicate(b,&sol);
611     
612     KrylovMinimize(A, b, x2);
613
614
615
616     MatMult(A,x2,sol);
617     VecAXPY(sol,-1,b);
618     VecNorm(sol, NORM_2, &norm);
619     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Error Krylov Minimization %g\n",norm);
620
621   }
622
623 {
624
625     Vec x2;
626     Vec sol;
627     VecDuplicate(b,&x2);
628     VecDuplicate(b,&sol);
629     
630     KrylovMinimizeLSQR(A, b, x2);
631
632
633
634     MatMult(A,x2,sol);
635     VecAXPY(sol,-1,b);
636     VecNorm(sol, NORM_2, &norm);
637     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Error KrylovLSQR Minimization %g\n",norm);
638
639   }
640
641
642   /*
643      Free work space.  All PETSc objects should be destroyed when they
644      are no longer needed.
645   */
646   ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = VecDestroy(&b);CHKERRQ(ierr);
647   ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr);
648   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
649   PetscPreLoadEnd();
650   /* -----------------------------------------------------------
651                       End of linear solver loop
652      ----------------------------------------------------------- */
653
654   ierr = PetscFinalize();
655   return 0;
656 }
657