Actual source code: ex49.c


  2: static char help[] = "Tests SeqSBAIJ factorizations for different block sizes\n\n";

  4: #include <petscksp.h>

  6: int main(int argc,char **args)
  7: {
  8:   Vec            x,b,u;
  9:   Mat            A,A2;
 10:   KSP            ksp;
 11:   PetscRandom    rctx;
 12:   PetscReal      norm;
 13:   PetscInt       i,j,k,l,n = 27,its,bs = 2,Ii,J;
 14:   PetscBool      test_hermitian = PETSC_FALSE, convert = PETSC_FALSE;
 15:   PetscScalar    v;

 17:   PetscInitialize(&argc,&args,(char*)0,help);
 18:   PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
 19:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 20:   PetscOptionsGetBool(NULL,NULL,"-herm",&test_hermitian,NULL);
 21:   PetscOptionsGetBool(NULL,NULL,"-conv",&convert,NULL);

 23:   MatCreate(PETSC_COMM_SELF,&A);
 24:   MatSetSizes(A,n*bs,n*bs,PETSC_DETERMINE,PETSC_DETERMINE);
 25:   MatSetBlockSize(A,bs);
 26:   MatSetType(A,MATSEQSBAIJ);
 27:   MatSetFromOptions(A);
 28:   MatSeqSBAIJSetPreallocation(A,bs,n,NULL);
 29:   MatSeqBAIJSetPreallocation(A,bs,n,NULL);
 30:   MatSeqAIJSetPreallocation(A,n*bs,NULL);
 31:   MatMPIAIJSetPreallocation(A,n*bs,NULL,n*bs,NULL);

 33:   PetscRandomCreate(PETSC_COMM_SELF,&rctx);
 34:   for (i=0; i<n; i++) {
 35:     for (j=i; j<n; j++) {
 36:       PetscRandomGetValue(rctx,&v);
 37:       if (PetscRealPart(v) < .1 || i == j) {
 38:         for (k=0; k<bs; k++) {
 39:           for (l=0; l<bs; l++) {
 40:             Ii = i*bs + k;
 41:             J = j*bs + l;
 42:             PetscRandomGetValue(rctx,&v);
 43:             if (Ii == J) v = PetscRealPart(v+3*n*bs);
 44:             MatSetValue(A,Ii,J,v,INSERT_VALUES);
 45:             if (test_hermitian) {
 46:               MatSetValue(A,J,Ii,PetscConj(v),INSERT_VALUES);
 47:             } else {
 48:               MatSetValue(A,J,Ii,v,INSERT_VALUES);
 49:             }
 50:           }
 51:         }
 52:       }
 53:     }
 54:   }
 55:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 56:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 58:   /* With complex numbers:
 59:      - PETSc cholesky does not support hermitian matrices
 60:      - CHOLMOD only supports hermitian matrices
 61:      - SUPERLU_DIST seems supporting both
 62:   */
 63:   if (test_hermitian) {
 64:     MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
 65:   }

 67:   {
 68:     Mat M;
 69:     MatComputeOperator(A,MATAIJ,&M);
 70:     MatViewFromOptions(M,NULL,"-expl_view");
 71:     MatDestroy(&M);
 72:   }

 74:   A2 = NULL;
 75:   if (convert) {
 76:     MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&A2);
 77:   }

 79:   VecCreate(PETSC_COMM_SELF,&u);
 80:   VecSetSizes(u,PETSC_DECIDE,n*bs);
 81:   VecSetFromOptions(u);
 82:   VecDuplicate(u,&b);
 83:   VecDuplicate(b,&x);

 85:   VecSet(u,1.0);
 86:   MatMult(A,u,b);

 88:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 89:                 Create the linear solver and set various options
 90:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 92:   /*
 93:      Create linear solver context
 94:   */
 95:   KSPCreate(PETSC_COMM_SELF,&ksp);

 97:   /*
 98:      Set operators.
 99:   */
100:   KSPSetOperators(ksp,A2 ? A2 : A,A);

102:   KSPSetFromOptions(ksp);

104:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105:                       Solve the linear system
106:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

108:   KSPSolve(ksp,b,x);

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:                       Check solution and clean up
112:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

114:   /*
115:      Check the error
116:   */
117:   VecAXPY(x,-1.0,u);
118:   VecNorm(x,NORM_2,&norm);
119:   KSPGetIterationNumber(ksp,&its);

121:   /*
122:      Print convergence information.  PetscPrintf() produces a single
123:      print statement from all processes that share a communicator.
124:      An alternative is PetscFPrintf(), which prints to a file.
125:   */
126:   if (norm > 100*PETSC_SMALL) {
127:     PetscPrintf(PETSC_COMM_SELF,"Norm of residual %g iterations %D bs %D\n",(double)norm,its,bs);
128:   }

130:   /*
131:      Free work space.  All PETSc objects should be destroyed when they
132:      are no longer needed.
133:   */
134:   KSPDestroy(&ksp);
135:   VecDestroy(&u);
136:   VecDestroy(&x);
137:   VecDestroy(&b);
138:   MatDestroy(&A);
139:   MatDestroy(&A2);
140:   PetscRandomDestroy(&rctx);

142:   /*
143:      Always call PetscFinalize() before exiting a program.  This routine
144:        - finalizes the PETSc libraries as well as MPI
145:        - provides summary and diagnostic information if certain runtime
146:          options are chosen (e.g., -log_view).
147:   */
148:   PetscFinalize();
149:   return 0;
150: }

152: /*TEST

154:    test:
155:       args: -mat_type {{aij baij sbaij}} -bs {{1 2 3 4 5 6 7 8 9 10 11 12}} -pc_type cholesky -herm 0 -conv {{0 1}}

157:    test:
158:       nsize: {{1 4}}
159:       suffix: cholmod
160:       requires: suitesparse
161:       args: -mat_type {{aij sbaij}} -bs 1 -pc_type cholesky -pc_factor_mat_solver_type cholmod -herm -conv {{0 1}}

163:    test:
164:       nsize: {{1 4}}
165:       suffix: superlu_dist
166:       requires: superlu_dist
167:       output_file: output/ex49_cholmod.out
168:       args: -mat_type mpiaij -bs 3 -pc_type cholesky -pc_factor_mat_solver_type superlu_dist -herm {{0 1}} -conv {{0 1}}

170:    test:
171:       suffix: mkl_pardiso
172:       requires: mkl_pardiso
173:       output_file: output/ex49_1.out
174:       args: -bs {{1 3}} -pc_type cholesky -pc_factor_mat_solver_type mkl_pardiso

176:    test:
177:       suffix: cg
178:       requires: complex
179:       output_file: output/ex49_cg.out
180:       args: -herm -ksp_cg_type hermitian -mat_type aij -ksp_type cg -pc_type jacobi -ksp_rtol 4e-07

182:    test:
183:       suffix: pipecg2
184:       requires: complex
185:       output_file: output/ex49_pipecg2.out
186:       args: -herm -mat_type aij -ksp_type pipecg2 -pc_type jacobi -ksp_rtol 4e-07 -ksp_norm_type {{preconditioned unpreconditioned natural}}

188: TEST*/