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;
 13:   PetscBool      flg = PETSC_FALSE;
 14:   PetscScalar    v;
 15:   PC             pc;
 16:   PetscInt       in;
 17:   Mat            F,B;
 18:   PetscBool      solve=PETSC_FALSE,sameA=PETSC_FALSE,setfromoptions_first=PETSC_FALSE;
 19: #if defined(PETSC_USE_LOG)
 20:   PetscLogStage stage;
 21: #endif
 22: #if !defined(PETSC_HAVE_MUMPS)
 23:   PetscMPIInt    size;
 24: #endif

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

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

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

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

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

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

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

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

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

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

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

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

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

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

159:   PetscFinalize();
160:   return 0;
161: }

163: /*TEST

165:     build:
166:       requires: !complex

168:     test:
169:       args:

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

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

179: TEST*/