Actual source code: ex17.c


  2: static char help[] = "Solves a linear system with KSP.  This problem is\n\
  3: intended to test the complex numbers version of various solvers.\n\n";

  5: #include <petscksp.h>

  7: typedef enum {TEST_1,TEST_2,TEST_3,HELMHOLTZ_1,HELMHOLTZ_2} TestType;
  8: extern PetscErrorCode FormTestMatrix(Mat,PetscInt,TestType);

 10: int main(int argc,char **args)
 11: {
 12:   Vec            x,b,u;      /* approx solution, RHS, exact solution */
 13:   Mat            A;            /* linear system matrix */
 14:   KSP            ksp;         /* KSP context */
 15:   PetscInt       n    = 10,its, dim,p = 1,use_random;
 16:   PetscScalar    none = -1.0,pfive = 0.5;
 17:   PetscReal      norm;
 18:   PetscRandom    rctx;
 19:   TestType       type;
 20:   PetscBool      flg;

 22:   PetscInitialize(&argc,&args,(char*)0,help);
 23:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 24:   PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
 25:   switch (p) {
 26:   case 1:  type = TEST_1;      dim = n;   break;
 27:   case 2:  type = TEST_2;      dim = n;   break;
 28:   case 3:  type = TEST_3;      dim = n;   break;
 29:   case 4:  type = HELMHOLTZ_1; dim = n*n; break;
 30:   case 5:  type = HELMHOLTZ_2; dim = n*n; break;
 31:   default: type = TEST_1;      dim = n;
 32:   }

 34:   /* Create vectors */
 35:   VecCreate(PETSC_COMM_WORLD,&x);
 36:   VecSetSizes(x,PETSC_DECIDE,dim);
 37:   VecSetFromOptions(x);
 38:   VecDuplicate(x,&b);
 39:   VecDuplicate(x,&u);

 41:   use_random = 1;
 42:   flg        = PETSC_FALSE;
 43:   PetscOptionsGetBool(NULL,NULL,"-norandom",&flg,NULL);
 44:   if (flg) {
 45:     use_random = 0;
 46:     VecSet(u,pfive);
 47:   } else {
 48:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
 49:     PetscRandomSetFromOptions(rctx);
 50:     VecSetRandom(u,rctx);
 51:   }

 53:   /* Create and assemble matrix */
 54:   MatCreate(PETSC_COMM_WORLD,&A);
 55:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
 56:   MatSetFromOptions(A);
 57:   MatSetUp(A);
 58:   FormTestMatrix(A,n,type);
 59:   MatMult(A,u,b);
 60:   flg  = PETSC_FALSE;
 61:   PetscOptionsGetBool(NULL,NULL,"-printout",&flg,NULL);
 62:   if (flg) {
 63:     MatView(A,PETSC_VIEWER_STDOUT_WORLD);
 64:     VecView(u,PETSC_VIEWER_STDOUT_WORLD);
 65:     VecView(b,PETSC_VIEWER_STDOUT_WORLD);
 66:   }

 68:   /* Create KSP context; set operators and options; solve linear system */
 69:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 70:   KSPSetOperators(ksp,A,A);
 71:   KSPSetFromOptions(ksp);
 72:   KSPSolve(ksp,b,x);
 73:   /* KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD); */

 75:   /* Check error */
 76:   VecAXPY(x,none,u);
 77:   VecNorm(x,NORM_2,&norm);
 78:   KSPGetIterationNumber(ksp,&its);
 79:   if (norm >= 1.e-12) {
 80:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
 81:   } else {
 82:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error < 1.e-12, Iterations %D\n",its);
 83:   }

 85:   /* Free work space */
 86:   VecDestroy(&x)); PetscCall(VecDestroy(&u);
 87:   VecDestroy(&b)); PetscCall(MatDestroy(&A);
 88:   if (use_random) PetscRandomDestroy(&rctx);
 89:   KSPDestroy(&ksp);
 90:   PetscFinalize();
 91:   return 0;
 92: }

 94: PetscErrorCode FormTestMatrix(Mat A,PetscInt n,TestType type)
 95: {
 96:   PetscScalar val[5];
 97:   PetscInt    i,j,Ii,J,col[5],Istart,Iend;

 99:   MatGetOwnershipRange(A,&Istart,&Iend);
100:   if (type == TEST_1) {
101:     val[0] = 1.0; val[1] = 4.0; val[2] = -2.0;
102:     for (i=1; i<n-1; i++) {
103:       col[0] = i-1; col[1] = i; col[2] = i+1;
104:       MatSetValues(A,1,&i,3,col,val,INSERT_VALUES);
105:     }
106:     i    = n-1; col[0] = n-2; col[1] = n-1;
107:     MatSetValues(A,1,&i,2,col,val,INSERT_VALUES);
108:     i    = 0; col[0] = 0; col[1] = 1; val[0] = 4.0; val[1] = -2.0;
109:     MatSetValues(A,1,&i,2,col,val,INSERT_VALUES);
110:   } else if (type == TEST_2) {
111:     val[0] = 1.0; val[1] = 0.0; val[2] = 2.0; val[3] = 1.0;
112:     for (i=2; i<n-1; i++) {
113:       col[0] = i-2; col[1] = i-1; col[2] = i; col[3] = i+1;
114:       MatSetValues(A,1,&i,4,col,val,INSERT_VALUES);
115:     }
116:     i    = n-1; col[0] = n-3; col[1] = n-2; col[2] = n-1;
117:     MatSetValues(A,1,&i,3,col,val,INSERT_VALUES);
118:     i    = 1; col[0] = 0; col[1] = 1; col[2] = 2;
119:     MatSetValues(A,1,&i,3,col,&val[1],INSERT_VALUES);
120:     i    = 0;
121:     MatSetValues(A,1,&i,2,col,&val[2],INSERT_VALUES);
122:   } else if (type == TEST_3) {
123:     val[0] = PETSC_i * 2.0;
124:     val[1] = 4.0; val[2] = 0.0; val[3] = 1.0; val[4] = 0.7;
125:     for (i=1; i<n-3; i++) {
126:       col[0] = i-1; col[1] = i; col[2] = i+1; col[3] = i+2; col[4] = i+3;
127:       MatSetValues(A,1,&i,5,col,val,INSERT_VALUES);
128:     }
129:     i    = n-3; col[0] = n-4; col[1] = n-3; col[2] = n-2; col[3] = n-1;
130:     MatSetValues(A,1,&i,4,col,val,INSERT_VALUES);
131:     i    = n-2; col[0] = n-3; col[1] = n-2; col[2] = n-1;
132:     MatSetValues(A,1,&i,3,col,val,INSERT_VALUES);
133:     i    = n-1; col[0] = n-2; col[1] = n-1;
134:     MatSetValues(A,1,&i,2,col,val,INSERT_VALUES);
135:     i    = 0; col[0] = 0; col[1] = 1; col[2] = 2; col[3] = 3;
136:     MatSetValues(A,1,&i,4,col,&val[1],INSERT_VALUES);
137:   } else if (type == HELMHOLTZ_1) {
138:     /* Problem domain: unit square: (0,1) x (0,1)
139:        Solve Helmholtz equation:
140:           -delta u - sigma1*u + i*sigma2*u = f,
141:            where delta = Laplace operator
142:        Dirichlet b.c.'s on all sides
143:      */
144:     PetscRandom rctx;
145:     PetscReal   h2,sigma1 = 5.0;
146:     PetscScalar sigma2;
147:     PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);
148:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
149:     PetscRandomSetFromOptions(rctx);
150:     PetscRandomSetInterval(rctx,0.0,PETSC_i);
151:     h2   = 1.0/((n+1)*(n+1));
152:     for (Ii=Istart; Ii<Iend; Ii++) {
153:       *val = -1.0; i = Ii/n; j = Ii - i*n;
154:       if (i>0) {
155:         J = Ii-n; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
156:       }
157:       if (i<n-1) {
158:         J = Ii+n; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
159:       }
160:       if (j>0) {
161:         J = Ii-1; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
162:       }
163:       if (j<n-1) {
164:         J = Ii+1; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
165:       }
166:       PetscRandomGetValue(rctx,&sigma2);
167:       *val = 4.0 - sigma1*h2 + sigma2*h2;
168:       MatSetValues(A,1,&Ii,1,&Ii,val,ADD_VALUES);
169:     }
170:     PetscRandomDestroy(&rctx);
171:   } else if (type == HELMHOLTZ_2) {
172:     /* Problem domain: unit square: (0,1) x (0,1)
173:        Solve Helmholtz equation:
174:           -delta u - sigma1*u = f,
175:            where delta = Laplace operator
176:        Dirichlet b.c.'s on 3 sides
177:        du/dn = i*alpha*u on (1,y), 0<y<1
178:      */
179:     PetscReal   h2,sigma1 = 200.0;
180:     PetscScalar alpha_h;
181:     PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);
182:     h2      = 1.0/((n+1)*(n+1));
183:     alpha_h = (PETSC_i * 10.0) / (PetscReal)(n+1);  /* alpha_h = alpha * h */
184:     for (Ii=Istart; Ii<Iend; Ii++) {
185:       *val = -1.0; i = Ii/n; j = Ii - i*n;
186:       if (i>0) {
187:         J = Ii-n; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
188:       }
189:       if (i<n-1) {
190:         J = Ii+n; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
191:       }
192:       if (j>0) {
193:         J = Ii-1; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
194:       }
195:       if (j<n-1) {
196:         J = Ii+1; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
197:       }
198:       *val = 4.0 - sigma1*h2;
199:       if (!((Ii+1)%n)) *val += alpha_h;
200:       MatSetValues(A,1,&Ii,1,&Ii,val,ADD_VALUES);
201:     }
202:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER_INPUT,"FormTestMatrix: unknown test matrix type");

204:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
205:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

207:   return 0;
208: }

210: /*TEST

212:     build:
213:       requires: complex

215:     test:
216:       args: -ksp_gmres_cgs_refinement_type refine_always -n 6 -ksp_monitor_short -p 5 -norandom -ksp_type gmres -pc_type jacobi -ksp_max_it 15
217:       requires: complex

219:     test:
220:       suffix: 2
221:       nsize: 3
222:       requires: complex
223:       args: -ksp_gmres_cgs_refinement_type refine_always -n 6 -ksp_monitor_short -p 5 -norandom -ksp_type gmres -pc_type jacobi -ksp_max_it 15
224:       output_file: output/ex17_1.out

226:     test:
227:       suffix: superlu_dist
228:       requires: superlu_dist complex
229:       args: -n 6 -p 5 -norandom -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_colperm MMD_ATA

231:     test:
232:       suffix: superlu_dist_2
233:       requires: superlu_dist complex
234:       nsize: 3
235:       output_file: output/ex17_superlu_dist.out
236:       args: -n 6 -p 5 -norandom -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_colperm MMD_ATA

238: TEST*/