Actual source code: ex8.c


  2: static char help[] = "Solves a linear system in parallel with KSP. \n\
  3: Contributed by Jose E. Roman, SLEPc developer, for testing repeated call of KSPSetOperators(), 2014 \n\n";

  5: #include <petscksp.h>
  6: int main(int argc,char **args)
  7: {
  8:   Vec            x,b,u;    /* approx solution, RHS, exact solution */
  9:   Mat            A;        /* linear system matrix */
 10:   KSP            ksp;      /* linear solver context */
 11:   PetscRandom    rctx;     /* random number generator context */
 12:   PetscInt       i,j,Ii,J,Istart,Iend,m = 8,n = 7;
 14:   PetscBool      flg = PETSC_FALSE;
 15:   PetscScalar    v;
 16:   PC             pc;
 17:   PetscInt       in;
 18:   Mat            F,B;
 19:   PetscBool      solve=PETSC_FALSE,sameA=PETSC_FALSE,setfromoptions_first=PETSC_FALSE;
 20: #if defined(PETSC_USE_LOG)
 21:   PetscLogStage stage;
 22: #endif
 23: #if !defined(PETSC_HAVE_MUMPS)
 24:   PetscMPIInt    size;
 25: #endif

 27:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 28:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 29:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 30:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 31:          Compute the matrix and right-hand-side vector that define
 32:          the linear system, Ax = b.
 33:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 34:   MatCreate(PETSC_COMM_WORLD,&A);
 35:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 36:   MatSetFromOptions(A);
 37:   MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
 38:   MatSeqAIJSetPreallocation(A,5,NULL);
 39:   MatSetUp(A);

 41:   MatGetOwnershipRange(A,&Istart,&Iend);

 43:   PetscLogStageRegister("Assembly", &stage);
 44:   PetscLogStagePush(stage);
 45:   for (Ii=Istart; Ii<Iend; Ii++) {
 46:     v = -1.0; i = Ii/n; j = Ii - i*n;
 47:     if (i>0)   {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 48:     if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 49:     if (j>0)   {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 50:     if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 51:     v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
 52:   }
 53:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 54:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 55:   PetscLogStagePop();

 57:   /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
 58:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);

 60:   /* Create parallel vectors. */
 61:   VecCreate(PETSC_COMM_WORLD,&u);
 62:   VecSetSizes(u,PETSC_DECIDE,m*n);
 63:   VecSetFromOptions(u);
 64:   VecDuplicate(u,&b);
 65:   VecDuplicate(b,&x);

 67:   /*
 68:      Set exact solution; then compute right-hand-side vector.
 69:      By default we use an exact solution of a vector with all
 70:      elements of 1.0;  Alternatively, using the runtime option
 71:      -random_sol forms a solution vector with random components.
 72:   */
 73:   PetscOptionsGetBool(NULL,NULL,"-random_exact_sol",&flg,NULL);
 74:   if (flg) {
 75:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
 76:     PetscRandomSetFromOptions(rctx);
 77:     VecSetRandom(u,rctx);
 78:     PetscRandomDestroy(&rctx);
 79:   } else {
 80:     VecSet(u,1.0);
 81:   }
 82:   MatMult(A,u,b);

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85:                 Create the linear solver and set various options
 86:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 87:   /* Create linear solver context */
 88:   KSPCreate(PETSC_COMM_WORLD,&ksp);

 90:   /* Set operators. */
 91:   KSPSetOperators(ksp,A,A);

 93:   KSPSetTolerances(ksp,1.e-2/((m+1)*(n+1)),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);

 95:   PetscOptionsGetBool(NULL,NULL,"-setfromoptions_first",&setfromoptions_first,NULL);
 96:   if (setfromoptions_first) {
 97:     /* code path for changing from KSPLSQR to KSPREONLY */
 98:     KSPSetFromOptions(ksp);
 99:   }
100:   KSPSetType(ksp,KSPPREONLY);
101:   KSPGetPC(ksp, &pc);
102:   PCSetType(pc,PCCHOLESKY);
103: #if defined(PETSC_HAVE_MUMPS)
104: #if defined(PETSC_USE_COMPLEX)
105:   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Spectrum slicing with MUMPS is not available for complex scalars");
106: #endif
107:   PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
108:   /*
109:      must use runtime option '-mat_mumps_icntl_13 1' (turn off ScaLAPACK for
110:      matrix inertia), currently there is no better way of setting this in program
111:   */
112:   PetscOptionsInsertString(NULL,"-mat_mumps_icntl_13 1");
113: #else
114:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
115:   if (size>1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Configure with MUMPS if you want to run this example in parallel");
116: #endif

118:   if (!setfromoptions_first) {
119:     /* when -setfromoptions_first is true, do not call KSPSetFromOptions() again and stick to KSPPREONLY */
120:     KSPSetFromOptions(ksp);
121:   }

123:   /* get inertia */
124:   PetscOptionsGetBool(NULL,NULL,"-solve",&solve,NULL);
125:   PetscOptionsGetBool(NULL,NULL,"-sameA",&sameA,NULL);
126:   KSPSetUp(ksp);
127:   PCFactorGetMatrix(pc,&F);
128:   MatGetInertia(F,&in,NULL,NULL);
129:   PetscPrintf(PETSC_COMM_WORLD,"INERTIA=%D\n",in);
130:   if (solve) {
131:     PetscPrintf(PETSC_COMM_WORLD,"Solving the intermediate KSP\n");
132:     KSPSolve(ksp,b,x);
133:   } else {PetscPrintf(PETSC_COMM_WORLD,"NOT Solving the intermediate KSP\n");}

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:                       Solve the linear system
137:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138:   MatDuplicate(A,MAT_COPY_VALUES,&B);
139:   if (sameA) {
140:     PetscPrintf(PETSC_COMM_WORLD,"Setting A\n");
141:     MatAXPY(A,1.1,B,DIFFERENT_NONZERO_PATTERN);
142:     KSPSetOperators(ksp,A,A);
143:   } else {
144:     PetscPrintf(PETSC_COMM_WORLD,"Setting B\n");
145:     MatAXPY(B,1.1,A,DIFFERENT_NONZERO_PATTERN);
146:     KSPSetOperators(ksp,B,B);
147:   }
148:   KSPSetUp(ksp);
149:   PCFactorGetMatrix(pc,&F);
150:   MatGetInertia(F,&in,NULL,NULL);
151:   PetscPrintf(PETSC_COMM_WORLD,"INERTIA=%D\n",in);
152:   KSPSolve(ksp,b,x);
153:   MatDestroy(&B);

155:   /* Free work space.*/
156:   KSPDestroy(&ksp);
157:   VecDestroy(&u);  VecDestroy(&x);
158:   VecDestroy(&b);  MatDestroy(&A);

160:   PetscFinalize();
161:   return ierr;
162: }

164: /*TEST

166:     build:
167:       requires: !complex

169:     test:
170:       args:

172:     test:
173:       suffix: 2
174:       args: -sameA

176:     test:
177:       suffix: 3
178:       args: -ksp_lsqr_monitor -ksp_type lsqr -setfromoptions_first {{0 1}separate output}

180: TEST*/