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;
 17:   PetscTruth     flg;
 18:   IS             perm,iperm;
 19:   MatFactorInfo  info;
 20:   PetscRandom    rand;

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

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

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

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

 57:   /* Test MatDuplicate() */
 58:   MatDuplicate(C,MAT_COPY_VALUES,&C1);

 60:   /* Test LU Factorization */
 61:   MatGetOrdering(C1,MATORDERING_NATURAL,&perm,&iperm);
 62:   MatLUFactorSymbolic(C1,perm,iperm,&info,&F);
 63:   for (nfact = 0; nfact < 2; nfact++){
 64:     if (!rank) printf(" LU nfact %d\n",nfact);
 65:     MatLUFactorNumeric(C1,&info,&F);

 67:     /* Test MatSolve() */
 68:     for (nsolve = 0; nsolve < 5; nsolve++){
 69:       VecGetArray(x,&array);
 70:       for (i=0; i<m; i++){
 71:         PetscRandomGetValue(rand,&rval);
 72:         array[i] = rval;
 73:       }
 74:       VecRestoreArray(x,&array);
 75:       VecCopy(x,u);
 76:       MatMult(C,x,b);

 78:       MatSolve(F,b,x);

 80:       /* Check the error */
 81:       VecAXPY(u,-1.0,x);  /* u <- (-1.0)x + u */
 82:       VecNorm(u,NORM_2,&norm);
 83:       if (!rank){
 84:         PetscPrintf(PETSC_COMM_SELF,"Norm of error %A\n",norm);
 85:       }
 86:     }
 87:   }
 88:   MatDestroy(C1);
 89:   MatDestroy(F);

 91:   /* Test Cholesky Factorization */
 92:   MatTranspose(C,&C1); /* C1 = C^T */
 93:   MatAXPY(C,1.0,C1,SAME_NONZERO_PATTERN); /* make C symmetric: C <- C + C^T */
 94:   MatShift(C,M);  /* make C positive definite */
 95:   MatDestroy(C1);
 96: 
 97:   MatSetOption(C,MAT_SYMMETRIC);
 98:   MatSetOption(C,MAT_SYMMETRY_ETERNAL);

100:   MatDuplicate(C,MAT_COPY_VALUES,&C1);
101:   MatCholeskyFactorSymbolic(C,perm,&info,&F);
102:   for (nfact = 0; nfact < 2; nfact++){
103:     if (!rank) printf(" Cholesky nfact %d\n",nfact);
104:     MatCholeskyFactorNumeric(C1,&info,&F);

106:     /* Test MatSolve() */
107:     for (nsolve = 0; nsolve < 5; nsolve++){
108:       VecGetArray(x,&array);
109:       for (i=0; i<m; i++){
110:         PetscRandomGetValue(rand,&rval);
111:         array[i] = rval;
112:       }
113:       VecRestoreArray(x,&array);
114:       VecCopy(x,u);
115:       MatMult(C,x,b);

117:       MatSolve(F,b,x);

119:       /* Check the error */
120:       VecAXPY(u,-1.0,x);  /* u <- (-1.0)x + u */
121:       VecNorm(u,NORM_2,&norm);
122:       if (!rank){
123:         PetscPrintf(PETSC_COMM_SELF,"Norm of error %A\n",norm);
124:       }
125:     }
126:   }
127:   MatDestroy(C1);
128:   MatDestroy(F);

130:   /* Free data structures */
131:   PetscRandomDestroy(rand);
132:   ISDestroy(perm);
133:   ISDestroy(iperm);
134:   VecDestroy(x);
135:   VecDestroy(b);
136:   VecDestroy(u);
137:   MatDestroy(C);

139: #else
140:   if (!rank) printf("This example needs PLAPLAPACK\n");
141: #endif
142:   PetscFinalize();
143:   return 0;
144: }