PETSc_DKE_ft

C/C++ source

0001 
0002 static char help[] = "Reads a PETSc matrix and vector from a binary file and solves a linear system.Saves a binary file.\n\
0003   See the 'Performance Hints' chapter of the\n\
0004 users manual for a discussion of preloading.  Input parameters include\n\
0005   -f0  : first file to load (small system)\n\
0006   -f1  : second file to load (larger system)\n\n\
0007   -trans  : solve transpose system instead\n\n";
0008 /*
0009   This code can be used to test PETSc interface to other packages.\n\
0010   Examples of command line options:       \n\
0011    ex10 -f0 <datafile> -ksp_type preonly  \n\
0012         -help -ksp_view                  \n\
0013         -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
0014         -ksp_type preonly -pc_type lu -mat_type aijspooles/superlu/superlu_dist/aijmumps \n\
0015         -ksp_type preonly -pc_type cholesky -mat_type sbaijspooles/dscpack/sbaijmumps \n\
0016         -f0 <A> -fB <B> -mat_type sbaijmumps -ksp_type preonly -pc_type cholesky -test_inertia -mat_sigma <sigma> \n\
0017    mpirun -np <np> ex10 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
0018  \n\n";
0019 */
0020 /*T
0021    Concepts: KSP^solving a linear system
0022    Processors: n
0023 T*/
0024 
0025 /* 
0026   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
0027   automatically includes:
0028      petsc.h       - base PETSc routines   petscvec.h - vectors
0029      petscsys.h    - system routines       petscmat.h - matrices
0030      petscis.h     - index sets            petscksp.h - Krylov subspace methods
0031      petscviewer.h - viewers               petscpc.h  - preconditioners
0032 */
0033 #include "petscksp.h"
0034 
0035 #undef __FUNCT__
0036 #define __FUNCT__ "main"
0037 int main(int argc,char **args)
0038 {
0039   KSP            ksp;             /* linear solver context */
0040   Mat            A,B;            /* matrix */
0041   Vec            x,b,u;          /* approx solution, RHS, exact solution */
0042   PetscViewer    fd;               /* viewer */
0043   char           file[3][PETSC_MAX_PATH_LEN];    /* input file name */
0044   char           outfile[PETSC_MAX_PATH_LEN];
0045   PetscTruth     table,flg,flgB=PETSC_FALSE,trans=PETSC_FALSE,partition=PETSC_FALSE;
0046   PetscErrorCode ierr;
0047   PetscInt       its,num_numfac;
0048   PetscReal      norm;
0049   PetscLogDouble tsetup,tsetup1,tsetup2,tsolve,tsolve1,tsolve2;
0050   PetscTruth     preload=PETSC_TRUE,diagonalscale,isSymmetric,cknorm=PETSC_FALSE,Test_MatDuplicate=PETSC_FALSE;
0051   PetscMPIInt    rank;
0052   PetscScalar    sigma;
0053   int size;
0054   
0055   PetscInitialize(&argc,&args,(char *)0,help);
0056   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
0057   ierr = PetscOptionsHasName(PETSC_NULL,"-table",&table);CHKERRQ(ierr);
0058   ierr = PetscOptionsHasName(PETSC_NULL,"-trans",&trans);CHKERRQ(ierr);
0059   ierr = PetscOptionsHasName(PETSC_NULL,"-partition",&partition);CHKERRQ(ierr);
0060  
0061 
0062   /* 
0063      Determine files from which we read the two linear systems
0064      (matrix and right-hand-side vector).
0065   */
0066   ierr = PetscOptionsGetString(PETSC_NULL,"-f",file[0],PETSC_MAX_PATH_LEN-1,&flg);CHKERRQ(ierr);
0067   if (flg) {
0068     ierr = PetscStrcpy(file[1],file[0]);CHKERRQ(ierr);
0069     preload = PETSC_FALSE;
0070   } else {
0071      SETERRQ(1,"Must indicate binary file with the -f option");
0072   }
0073 
0074   /* -----------------------------------------------------------
0075                   Beginning of linear solver loop
0076      ----------------------------------------------------------- */
0077  
0078   PreLoadBegin(preload,"Load system");
0079 
0080     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
0081                            Load system
0082      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
0083 
0084     /* 
0085        Open binary file.  Note that we use FILE_MODE_READ to indicate
0086        reading from this file.
0087     */
0088     ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PreLoadIt],FILE_MODE_READ,&fd);CHKERRQ(ierr);
0089     
0090     /*
0091        Load the matrix and vector; then destroy the viewer.
0092     */
0093     ierr = MatLoad(fd,MATAIJ,&A);CHKERRQ(ierr);
0094     if (!preload){
0095       flg = PETSC_FALSE;
0096       ierr = PetscOptionsGetString(PETSC_NULL,"-rhs",file[2],PETSC_MAX_PATH_LEN-1,&flg);CHKERRQ(ierr);
0097       if (flg){ /* rhs is stored in a separate file */
0098         ierr = PetscViewerDestroy(fd);CHKERRQ(ierr); 
0099         ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);CHKERRQ(ierr);
0100       } else{SETERRQ(1,"Must indicate rhs binary file with the -rhs option");}
0101     }
0102     if (rank){
0103         ierr = PetscExceptionTry1(VecLoad(fd,PETSC_NULL,&b),PETSC_ERR_FILE_UNEXPECTED);
0104     } else {
0105       ierr = PetscExceptionTry1(VecLoad(fd,PETSC_NULL,&b),PETSC_ERR_FILE_READ); 
0106     }   
0107     if (PetscExceptionCaught(ierr,PETSC_ERR_FILE_UNEXPECTED) || PetscExceptionCaught(ierr,PETSC_ERR_FILE_READ)) { /* if file contains no RHS, then use a vector of all ones */
0108       PetscInt    m;
0109       PetscScalar one = 1.0;
0110       ierr = PetscInfo(0,"Using vector of ones for RHS\n");CHKERRQ(ierr);
0111       ierr = MatGetLocalSize(A,&m,PETSC_NULL);CHKERRQ(ierr);
0112       ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
0113       ierr = VecSetSizes(b,m,PETSC_DECIDE);CHKERRQ(ierr);
0114       ierr = VecSetFromOptions(b);CHKERRQ(ierr);
0115       ierr = VecSet(b,one);CHKERRQ(ierr);
0116     } else CHKERRQ(ierr); 
0117     ierr = PetscViewerDestroy(fd);CHKERRQ(ierr); 
0118 
0119     /* Test MatDuplicate() */
0120     if (Test_MatDuplicate){
0121       ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
0122       ierr = MatEqual(A,B,&flg);CHKERRQ(ierr);
0123       if (!flg){
0124         PetscPrintf(PETSC_COMM_WORLD,"  A != B \n");CHKERRQ(ierr);
0125       } 
0126       ierr = MatDestroy(B);CHKERRQ(ierr); 
0127     }
0128 
0129     /* Add a shift to A */
0130     ierr = PetscOptionsGetScalar(PETSC_NULL,"-mat_sigma",&sigma,&flg);CHKERRQ(ierr);
0131     if(flg) {
0132       ierr = PetscOptionsGetString(PETSC_NULL,"-fB",file[2],PETSC_MAX_PATH_LEN-1,&flgB);CHKERRQ(ierr);
0133       if (flgB){
0134         /* load B to get A = A + sigma*B */
0135         ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);CHKERRQ(ierr);
0136         ierr  = MatLoad(fd,MATAIJ,&B);CHKERRQ(ierr);
0137         ierr = PetscViewerDestroy(fd);CHKERRQ(ierr);
0138         ierr = MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); /* A <- sigma*B + A */  
0139       } else {
0140         ierr = MatShift(A,sigma);CHKERRQ(ierr); 
0141       }
0142     }
0143 
0144     /* Make A singular for testing zero-pivot of ilu factorization        */
0145     /* Example: ./ex10 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_nonzero */
0146     ierr = PetscOptionsHasName(PETSC_NULL, "-test_zeropivot", &flg);CHKERRQ(ierr);
0147     if (flg) {
0148       PetscInt          row,ncols;
0149       const PetscInt    *cols;
0150       const PetscScalar *vals;
0151       PetscTruth        flg1=PETSC_FALSE;
0152       PetscScalar       *zeros;
0153       row = 0;      
0154       ierr = MatGetRow(A,row,&ncols,&cols,&vals);CHKERRQ(ierr);     
0155       ierr = PetscMalloc(sizeof(PetscScalar)*(ncols+1),&zeros);
0156       ierr = PetscMemzero(zeros,(ncols+1)*sizeof(PetscScalar));CHKERRQ(ierr);
0157       ierr = PetscOptionsHasName(PETSC_NULL, "-set_row_zero", &flg1);CHKERRQ(ierr);
0158       if (flg1){ /* set entire row as zero */
0159         ierr = MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);CHKERRQ(ierr);
0160       } else { /* only set (row,row) entry as zero */
0161         ierr = MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);CHKERRQ(ierr);
0162       }
0163       ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
0164       ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
0165     }
0166 
0167     /* Check whether A is symmetric */
0168     ierr = PetscOptionsHasName(PETSC_NULL, "-check_symmetry", &flg);CHKERRQ(ierr);
0169     if (flg) {
0170       Mat Atrans;
0171       ierr = MatTranspose(A, &Atrans);
0172       ierr = MatEqual(A, Atrans, &isSymmetric);
0173       if (isSymmetric) {
0174         PetscPrintf(PETSC_COMM_WORLD,"A is symmetric \n");CHKERRQ(ierr);
0175       } else {
0176         PetscPrintf(PETSC_COMM_WORLD,"A is non-symmetric \n");CHKERRQ(ierr);
0177       }
0178       ierr = MatDestroy(Atrans);CHKERRQ(ierr);
0179     }
0180 
0181     /* 
0182        If the loaded matrix is larger than the vector (due to being padded 
0183        to match the block size of the system), then create a new padded vector.
0184     */
0185     { 
0186       PetscInt    m,n,j,mvec,start,end,indx;
0187       Vec         tmp;
0188       PetscScalar *bold;
0189 
0190       /* Create a new vector b by padding the old one */
0191       ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
0192       if (m != n) {
0193         SETERRQ2(PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
0194       }
0195       ierr = VecCreate(PETSC_COMM_WORLD,&tmp);CHKERRQ(ierr);
0196       ierr = VecSetSizes(tmp,m,PETSC_DECIDE);CHKERRQ(ierr);
0197       ierr = VecSetFromOptions(tmp);CHKERRQ(ierr);
0198       ierr = VecGetOwnershipRange(b,&start,&end);CHKERRQ(ierr);
0199       ierr = VecGetLocalSize(b,&mvec);CHKERRQ(ierr);
0200       ierr = VecGetArray(b,&bold);CHKERRQ(ierr);
0201       for (j=0; j/* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
0216                       Setup solve for system
0217      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
0218 
0219 
0220     if (partition) {
0221       MatPartitioning mpart;
0222       IS              mis,nis,isn,is;
0223       PetscInt        *count;
0224       PetscMPIInt     size;
0225       Mat             BB;
0226       ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
0227       ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
0228       ierr = PetscMalloc(size*sizeof(PetscInt),&count);CHKERRQ(ierr);
0229       ierr = MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);CHKERRQ(ierr);
0230       ierr = MatPartitioningSetAdjacency(mpart, A);CHKERRQ(ierr);
0231       /* ierr = MatPartitioningSetVertexWeights(mpart, weight);CHKERRQ(ierr); */
0232       ierr = MatPartitioningSetFromOptions(mpart);CHKERRQ(ierr);
0233       ierr = MatPartitioningApply(mpart, &mis);CHKERRQ(ierr);
0234       ierr = MatPartitioningDestroy(mpart);CHKERRQ(ierr);
0235       ierr = ISPartitioningToNumbering(mis,&nis);CHKERRQ(ierr);
0236       ierr = ISPartitioningCount(mis,count);CHKERRQ(ierr);
0237       ierr = ISDestroy(mis);CHKERRQ(ierr);
0238       ierr = ISInvertPermutation(nis, count[rank], &is);CHKERRQ(ierr);
0239       ierr = PetscFree(count);CHKERRQ(ierr);
0240       ierr = ISDestroy(nis);CHKERRQ(ierr);
0241       ierr = ISSort(is);CHKERRQ(ierr);
0242       ierr = ISAllGather(is,&isn);CHKERRQ(ierr);
0243       ierr = MatGetSubMatrix(A,is,isn,PETSC_DECIDE,MAT_INITIAL_MATRIX,&BB);CHKERRQ(ierr);
0244 
0245       /* need to move the vector also */
0246       ierr = ISDestroy(is);CHKERRQ(ierr);
0247       ierr = ISDestroy(isn);CHKERRQ(ierr);
0248       ierr = MatDestroy(A);CHKERRQ(ierr);
0249       A    = BB;
0250     }
0251  
0252     /*
0253        Conclude profiling last stage; begin profiling next stage.
0254     */
0255     PreLoadStage("KSPSetUp");
0256 
0257     /*
0258        We also explicitly time this stage via PetscGetTime()
0259     */
0260     ierr = PetscGetTime(&tsetup1);CHKERRQ(ierr);
0261 
0262     /*
0263        Create linear solver; set operators; set runtime options.
0264     */
0265     ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
0266 
0267     num_numfac = 1;
0268     ierr = PetscOptionsGetInt(PETSC_NULL,"-num_numfac",&num_numfac,PETSC_NULL);CHKERRQ(ierr);
0269     while ( num_numfac-- ){
0270       /* ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); */
0271     ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
0272     ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
0273 
0274     /* 
0275        Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
0276        enable more precise profiling of setting up the preconditioner.
0277        These calls are optional, since both will be called within
0278        KSPSolve() if they haven't been called already.
0279     */
0280     ierr = KSPSetUp(ksp);CHKERRQ(ierr);
0281     ierr = KSPSetUpOnBlocks(ksp);CHKERRQ(ierr);
0282     ierr = PetscGetTime(&tsetup2);CHKERRQ(ierr);
0283     tsetup = tsetup2 - tsetup1;
0284 
0285     /*
0286       Test MatGetInertia()
0287       Usage:
0288       ex10 -f0 <mat_binaryfile> -ksp_type preonly -pc_type cholesky -mat_type seqsbaij -test_inertia -mat_sigma <sigma>
0289      */
0290     ierr = PetscOptionsHasName(PETSC_NULL,"-test_inertia",&flg);CHKERRQ(ierr);
0291     if (flg){
0292       PC        pc;
0293       PetscInt  nneg, nzero, npos;
0294       Mat       F;
0295       
0296       ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
0297       ierr = PCGetFactoredMatrix(pc,&F);CHKERRQ(ierr);
0298       ierr = MatGetInertia(F,&nneg,&nzero,&npos);CHKERRQ(ierr);
0299       ierr = PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
0300     }
0301 
0302     /*
0303        Tests "diagonal-scaling of preconditioned residual norm" as used 
0304        by many ODE integrator codes including SUNDIALS. Note this is different
0305        than diagonally scaling the matrix before computing the preconditioner
0306     */
0307     ierr = PetscOptionsHasName(PETSC_NULL,"-diagonal_scale",&diagonalscale);CHKERRQ(ierr);
0308     if (diagonalscale) {
0309       PC       pc;
0310       PetscInt j,start,end,n;
0311       Vec      scale;
0312       
0313       ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
0314       ierr = VecGetSize(x,&n);CHKERRQ(ierr);
0315       ierr = VecDuplicate(x,&scale);CHKERRQ(ierr);
0316       ierr = VecGetOwnershipRange(scale,&start,&end);CHKERRQ(ierr);
0317       for (j=start; j<end; j++) {
0318         ierr = VecSetValue(scale,j,((PetscReal)(j+1))/((PetscReal)n),INSERT_VALUES);CHKERRQ(ierr);
0319       }
0320       ierr = VecAssemblyBegin(scale);CHKERRQ(ierr);
0321       ierr = VecAssemblyEnd(scale);CHKERRQ(ierr);
0322       ierr = PCDiagonalScaleSet(pc,scale);CHKERRQ(ierr);
0323       ierr = VecDestroy(scale);CHKERRQ(ierr);
0324 
0325     }
0326 
0327     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
0328                            Solve system
0329      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
0330 
0331     /*
0332        Begin profiling next stage
0333     */
0334     PreLoadStage("KSPSolve");
0335 
0336     /*
0337        Solve linear system; we also explicitly time this stage.
0338     */
0339     ierr = PetscGetTime(&tsolve1);CHKERRQ(ierr);
0340     if (trans) {
0341       ierr = KSPSolveTranspose(ksp,b,x);CHKERRQ(ierr);
0342       ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
0343     } else {
0344       PetscInt  num_rhs=1;
0345       ierr = PetscOptionsGetInt(PETSC_NULL,"-num_rhs",&num_rhs,PETSC_NULL);CHKERRQ(ierr);
0346       ierr = PetscOptionsHasName(PETSC_NULL,"-cknorm",&cknorm);CHKERRQ(ierr);
0347       while ( num_rhs-- ) {
0348         ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
0349       }
0350       ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
0351       if (cknorm){   /* Check error for each rhs */
0352         if (trans) {
0353           ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr);
0354         } else {
0355           ierr = MatMult(A,x,u);CHKERRQ(ierr);
0356         }
0357         ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr);
0358         ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
0359         ierr = PetscPrintf(PETSC_COMM_WORLD,"  Number of iterations = %3D\n",its);CHKERRQ(ierr);
0360         ierr = PetscPrintf(PETSC_COMM_WORLD,"  Residual norm %A\n",norm);CHKERRQ(ierr);
0361        }
0362      
0363      
0364     } /* while ( num_rhs-- ) */
0365     ierr = PetscGetTime(&tsolve2);CHKERRQ(ierr);
0366     tsolve = tsolve2 - tsolve1;
0367 
0368    /* 
0369        Conclude profiling this stage
0370     */
0371     PreLoadStage("Cleanup");
0372 
0373     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
0374             Check error, print output, free data structures.
0375      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
0376 
0377     /* 
0378        Check error
0379     */
0380     if (trans) {
0381       ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr);
0382     } else {
0383       ierr = MatMult(A,x,u);CHKERRQ(ierr);
0384     }
0385     ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr);
0386     ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
0387 
0388     /*
0389        Write output (optinally using table for solver details).
0390         - PetscPrintf() handles output for multiprocessor jobs 
0391           by printing from only one processor in the communicator.
0392         - KSPView() prints information about the linear solver.
0393     */
0394     if (table) {
0395       char        *matrixname,kspinfo[120];
0396       PetscViewer viewer;
0397 
0398       /*
0399          Open a string viewer; then write info to it.
0400       */
0401       ierr = PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,120,&viewer);CHKERRQ(ierr);
0402       ierr = KSPView(ksp,viewer);CHKERRQ(ierr);
0403       ierr = PetscStrrchr(file[PreLoadIt],'/',&matrixname);CHKERRQ(ierr);
0404       ierr = PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %2.1e %2.1e %2.1e %s \n",
0405                 matrixname,its,norm,tsetup+tsolve,tsetup,tsolve,kspinfo);CHKERRQ(ierr);
0406       
0407       
0408      
0409       
0410       
0411       /*
0412          Destroy the viewer
0413       */
0414       ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
0415     } else {
0416       ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);CHKERRQ(ierr);
0417       ierr = PetscPrintf(PETSC_COMM_WORLD,"Residual norm %A\n",norm);CHKERRQ(ierr);
0418     }
0419 
0420     ierr = PetscOptionsHasName(PETSC_NULL, "-ksp_reason", &flg);CHKERRQ(ierr);
0421     if (flg){
0422       KSPConvergedReason reason;
0423       ierr = KSPGetConvergedReason(ksp,&reason);CHKERRQ(ierr);
0424       PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason); 
0425     }
0426        
0427     } /* while ( num_numfac-- ) */
0428 
0429 /*EDIT=================================================*/
0430 
0431 ierr = PetscOptionsGetString(PETSC_NULL,"-scratchdir",outfile,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
0432   if (flg) {
0433     ierr = PetscStrcat(outfile,"/PETSc_out_X.bin");CHKERRQ(ierr);
0434   } else { SETERRQ(1,"Must indicate scratch directory with the -scracthdir option")};
0435 
0436 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,outfile,FILE_MODE_WRITE,&fd);CHKERRQ(ierr);
0437 VecView(x,fd);
0438 PetscViewerDestroy(fd);
0439 
0440 /*EDIT=================================================*/
0441 
0442     /* 
0443        Free work space.  All PETSc objects should be destroyed when they
0444        are no longer needed.
0445     */
0446     ierr = VecGetSize(x,&size);CHKERRQ(ierr);
0447     printf("size=%i\n",size);
0448     ierr = VecGetSize(b,&size);CHKERRQ(ierr);
0449     printf("size=%i\n",size);
0450     ierr = MatDestroy(A);CHKERRQ(ierr); ierr = VecDestroy(b);CHKERRQ(ierr);
0451     ierr = VecDestroy(u);CHKERRQ(ierr); ierr = VecDestroy(x);CHKERRQ(ierr);
0452     ierr = KSPDestroy(ksp);CHKERRQ(ierr); 
0453     if (flgB) { ierr = MatDestroy(B);CHKERRQ(ierr); }
0454   PreLoadEnd();
0455   /* -----------------------------------------------------------
0456                       End of linear solver loop
0457      ----------------------------------------------------------- */
0458 
0459   ierr = PetscFinalize();CHKERRQ(ierr);
0460   return 0;
0461 }
0462