Actual source code: ex18.c
petsc-3.4.2 2013-07-02

      2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
      3: Input arguments are:\n\
      4:   -f <input_file> : file to load.  For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\n";

      6: #include <petscmat.h>
      7: #include <petscksp.h>
      8: #include <petsctime.h>

     12: int main(int argc,char **args)
     13: {
     15:   PetscInt       its,m,n,mvec;
     16:   PetscLogDouble time1,time2,time;
     17:   PetscReal      norm;
     18:   Vec            x,b,u;
     19:   Mat            A;
     20:   KSP            ksp;
     21:   char           file[PETSC_MAX_PATH_LEN];
     22:   PetscViewer    fd;
     23:   PetscLogStage  stage1;

     25:   PetscInitialize(&argc,&args,(char*)0,help);

     27:   /* Read matrix and RHS */
     28:   PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,NULL);
     29:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
     30:   MatCreate(PETSC_COMM_WORLD,&A);
     31:   MatSetType(A,MATSEQAIJ);
     32:   MatLoad(A,fd);
     33:   VecCreate(PETSC_COMM_WORLD,&b);
     34:   VecLoad(b,fd);
     35:   PetscViewerDestroy(&fd);

     37:   /*
     38:      If the load matrix is larger then the vector, due to being padded
     39:      to match the blocksize then create a new padded vector
     40:   */
     41:   MatGetSize(A,&m,&n);
     42:   VecGetSize(b,&mvec);
     43:   if (m > mvec) {
     44:     Vec         tmp;
     45:     PetscScalar *bold,*bnew;
     46:     /* create a new vector b by padding the old one */
     47:     VecCreate(PETSC_COMM_WORLD,&tmp);
     48:     VecSetSizes(tmp,PETSC_DECIDE,m);
     49:     VecSetFromOptions(tmp);
     50:     VecGetArray(tmp,&bnew);
     51:     VecGetArray(b,&bold);
     52:     PetscMemcpy(bnew,bold,mvec*sizeof(PetscScalar));
     53:     VecDestroy(&b);
     54:     b    = tmp;
     55:   }

     57:   /* Set up solution */
     58:   VecDuplicate(b,&x);
     59:   VecDuplicate(b,&u);
     60:   VecSet(x,0.0);

     62:   /* Solve system */
     63:   PetscLogStageRegister("Stage 1",&stage1);
     64:   PetscLogStagePush(stage1);
     65:   KSPCreate(PETSC_COMM_WORLD,&ksp);
     66:   KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
     67:   KSPSetFromOptions(ksp);
     68:   PetscTime(&time1);
     69:   KSPSolve(ksp,b,x);
     70:   PetscTime(&time2);
     71:   time = time2 - time1;
     72:   PetscLogStagePop();

     74:   /* Show result */
     75:   MatMult(A,x,u);
     76:   VecAXPY(u,-1.0,b);
     77:   VecNorm(u,NORM_2,&norm);
     78:   KSPGetIterationNumber(ksp,&its);
     79:   PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
     80:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
     81:   PetscPrintf(PETSC_COMM_WORLD,"Time for solve = %5.2f seconds\n",time);

     83:   /* Cleanup */
     84:   KSPDestroy(&ksp);
     85:   VecDestroy(&x);
     86:   VecDestroy(&b);
     87:   VecDestroy(&u);
     88:   MatDestroy(&A);

     90:   PetscFinalize();
     91:   return 0;
     92: }