Actual source code: ex103.c

  1: 
  2: static char help[] = "Tests PLAPACK interface.\n\n";

 4:  #include petscmat.h

  8: int main(int argc,char **args)
  9: {
 10:   Mat            C,C1,F;
 11:   Vec            u,x,b;
 13:   PetscMPIInt    rank,nproc;
 14:   PetscInt       i,M = 10,m,n,nfact,nsolve;
 15:   PetscScalar    *array,rval;
 16:   PetscReal      norm,tol=1.e-12;
 17:   IS             perm,iperm;
 18:   MatFactorInfo  info;
 19:   PetscRandom    rand;
 20:   PetscTruth     flg;

 22:   PetscInitialize(&argc,&args,(char *)0,help);
 23:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 24:   MPI_Comm_size(PETSC_COMM_WORLD, &nproc);

 26:   /* Create matrix and vectors */
 27:   PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
 28:   MatCreate(PETSC_COMM_WORLD,&C);
 29:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,M,M);
 30:   MatSetType(C,MATDENSE);
 31:   MatSetFromOptions(C);
 32: 
 33:   MatGetLocalSize(C,&m,&n);
 34:   if (m != n) SETERRQ2(PETSC_ERR_ARG_WRONG,"Matrix local size m %d must equal n %d",m,n);

 36:   VecCreate(PETSC_COMM_WORLD,&x);
 37:   VecSetSizes(x,n,PETSC_DECIDE);
 38:   VecSetFromOptions(x);
 39:   VecDuplicate(x,&b);
 40:   VecDuplicate(x,&u); /* save the true solution */

 42:   /* Assembly */
 43:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 44:   PetscRandomSetFromOptions(rand);
 45:   MatGetArray(C,&array);
 46:   for (i=0; i<m*M; i++){
 47:     PetscRandomGetValue(rand,&rval);
 48:     array[i] = rval;
 49:   }
 50:   MatRestoreArray(C,&array);
 51:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 52:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
 53:   /*if (!rank) {printf("main, C: \n");}
 54:     MatView(C,PETSC_VIEWER_STDOUT_WORLD); */

 56:   /* Test MatDuplicate() */
 57:   MatDuplicate(C,MAT_COPY_VALUES,&C1);
 58:   MatEqual(C,C1,&flg);
 59:   if (!flg){
 60:     SETERRQ(PETSC_ERR_ARG_WRONG,"Duplicate C1 != C");
 61:   }

 63:   /* Test LU Factorization */
 64:   MatGetOrdering(C1,MATORDERING_NATURAL,&perm,&iperm);
 65:   if (nproc == 1){
 66:     MatGetFactor(C1,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&F);
 67:   } else {
 68:     MatGetFactor(C1,MAT_SOLVER_PLAPACK,MAT_FACTOR_LU,&F);
 69:   }
 70:   MatLUFactorSymbolic(F,C1,perm,iperm,&info);

 72:   for (nfact = 0; nfact < 2; nfact++){
 73:     if (!rank) printf(" LU nfact %d\n",nfact);
 74:     MatLUFactorNumeric(F,C1,&info);

 76:     /* Test MatSolve() */
 77:     for (nsolve = 0; nsolve < 5; nsolve++){
 78:       VecGetArray(x,&array);
 79:       for (i=0; i<m; i++){
 80:         PetscRandomGetValue(rand,&rval);
 81:         array[i] = rval;
 82:       }
 83:       VecRestoreArray(x,&array);
 84:       VecCopy(x,u);
 85:       MatMult(C,x,b);

 87:       MatSolve(F,b,x);

 89:       /* Check the error */
 90:       VecAXPY(u,-1.0,x);  /* u <- (-1.0)x + u */
 91:       VecNorm(u,NORM_2,&norm);
 92:       if (norm > tol){
 93:         if (!rank){
 94:           PetscPrintf(PETSC_COMM_SELF,"Error: Norm of error %g, LU nfact %d\n",norm,nfact);
 95:         }
 96:       }
 97:     }
 98:   }
 99:   MatDestroy(C1);
100:   MatDestroy(F);

102:   /* Test Cholesky Factorization */
103:   MatTranspose(C,MAT_INITIAL_MATRIX,&C1); /* C1 = C^T */
104:   MatAXPY(C,1.0,C1,SAME_NONZERO_PATTERN); /* make C symmetric: C <- C + C^T */
105:   MatShift(C,M);  /* make C positive definite */
106:   MatDestroy(C1);
107: 
108:   MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);
109:   MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
110: 
111:   if (nproc == 1){
112:     MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&F);
113:   } else {
114:     MatGetFactor(C,MAT_SOLVER_PLAPACK,MAT_FACTOR_CHOLESKY,&F);
115:   }
116:   MatCholeskyFactorSymbolic(F,C,perm,&info);
117:   for (nfact = 0; nfact < 2; nfact++){
118:     if (!rank) printf(" Cholesky nfact %d\n",nfact);
119:     MatCholeskyFactorNumeric(F,C,&info);

121:     /* Test MatSolve() */
122:     for (nsolve = 0; nsolve < 5; nsolve++){
123:       VecGetArray(x,&array);
124:       for (i=0; i<m; i++){
125:         PetscRandomGetValue(rand,&rval);
126:         array[i] = rval;
127:       }
128:       VecRestoreArray(x,&array);
129:       VecCopy(x,u);
130:       MatMult(C,x,b);

132:       MatSolve(F,b,x);

134:       /* Check the error */
135:       VecAXPY(u,-1.0,x);  /* u <- (-1.0)x + u */
136:       VecNorm(u,NORM_2,&norm);
137:       if (norm > tol){
138:         if (!rank){
139:           PetscPrintf(PETSC_COMM_SELF,"Error: Norm of error %g, Cholesky nfact %d\n",norm,nfact);
140:         }
141:       }
142:     }
143:   }
144:   MatDestroy(F);

146:   /* Free data structures */
147:   PetscRandomDestroy(rand);
148:   ISDestroy(perm);
149:   ISDestroy(iperm);
150:   VecDestroy(x);
151:   VecDestroy(b);
152:   VecDestroy(u);
153:   MatDestroy(C);

155:   PetscFinalize();
156:   return 0;
157: }