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