Actual source code: ex17.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      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 */
     16:   PetscInt       n    = 10,its, dim,p = 1,use_random;
     17:   PetscScalar    none = -1.0,pfive = 0.5;
     18:   PetscReal      norm;
     19:   PetscRandom    rctx;
     20:   TestType       type;
     21:   PetscBool      flg;

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

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

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

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

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

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

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

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

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

    206:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    207:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

    209:   return 0;
    210: }

    212: /*TEST

    214:     build:
    215:       requires: complex

    217:     test:
    218:       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
    219:       requires: complex

    221:     test:
    222:       suffix: 2
    223:       nsize: 3
    224:       requires: complex
    225:       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
    226:       output_file: output/ex17_1.out

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

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

    240: TEST*/