Actual source code: ex103.c

petsc-3.3-p6 2013-02-11
  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,start,end;
 15:   PetscScalar    *array,rval;
 16:   PetscReal      norm,tol=1.e-12;
 17:   IS             perm,iperm;
 18:   MatFactorInfo  info;
 19:   PetscRandom    rand;
 20:   PetscBool      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 C */
 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:   MatSetUp(C);

 34:   /* Assembly C */
 35:   MatGetOwnershipRange(C,&start,&end);
 36:   m    = end - start;
 37:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 38:   PetscRandomSetFromOptions(rand);
 39:   MatGetArray(C,&array);
 40:   for (i=0; i<m*M; i++){
 41:     PetscRandomGetValue(rand,&rval);
 42:     array[i] = rval;
 43:   }
 44:   MatRestoreArray(C,&array);
 45:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 46:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
 47:   /* if (!rank) {printf("main, C: \n");}
 48:      MatView(C,PETSC_VIEWER_STDOUT_WORLD); */
 49: 
 50:   /* Create vectors */
 51:   MatGetLocalSize(C,&m,&n);
 52:   if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix local size m %d must equal n %d",m,n);
 53:   VecCreate(PETSC_COMM_WORLD,&x);
 54:   VecSetSizes(x,n,PETSC_DECIDE);
 55:   VecSetFromOptions(x);
 56:   VecDuplicate(x,&b);
 57:   VecDuplicate(x,&u); /* save the true solution */

 59:   /* Test MatDuplicate() */
 60:   MatDuplicate(C,MAT_COPY_VALUES,&C1);
 61:   MatEqual(C,C1,&flg);
 62:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Duplicate C1 != C");

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

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

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

 88:       MatSolve(F,b,x);

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

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

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

133:       MatSolve(F,b,x);

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

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

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