2: static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for a sequential dense matrix. \n\
  3:                       For MATSEQDENSE matrix, the factorization is just a thin wrapper to LAPACK \n\n";

  5: #include <petscmat.h>

  9: int main(int argc,char **argv)
 10: {
 11:   Mat            mat,F,RHS,SOLU;
 12:   MatInfo        info;
 14:   PetscInt       m = 10,n = 10,i,j,rstart,rend,nrhs=2;
 15:   PetscScalar    value = 1.0;
 16:   Vec            x,y,b,ytmp;
 17:   PetscReal      norm,tol=1.e-15;
 18:   PetscMPIInt    size;
 19:   PetscScalar    *rhs_array,*solu_array;
 20:   PetscRandom    rand;
 21:   PetscScalar    *array,rval;

 23:   PetscInitialize(&argc,&argv,(char*) 0,help);
 24:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 25:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");

 27:   /* create single vectors */
 28:   VecCreate(PETSC_COMM_WORLD,&y);
 29:   VecSetSizes(y,PETSC_DECIDE,m);
 30:   VecSetFromOptions(y);
 31:   VecDuplicate(y,&x);
 32:   VecDuplicate(y,&ytmp);
 33:   VecSet(x,value);
 34:   VecCreate(PETSC_COMM_WORLD,&b);
 35:   VecSetSizes(b,PETSC_DECIDE,n);
 36:   VecSetFromOptions(b);

 38:   /* create multiple vectors RHS and SOLU */
 39:   MatCreate(PETSC_COMM_WORLD,&RHS);
 40:   MatSetSizes(RHS,PETSC_DECIDE,PETSC_DECIDE,n,nrhs);
 41:   MatSetType(RHS,MATDENSE);
 42:   MatSetFromOptions(RHS);
 43:   MatSeqDenseSetPreallocation(RHS,NULL);

 45:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 46:   PetscRandomSetFromOptions(rand);
 47:   MatDenseGetArray(RHS,&array);
 48:   for (j=0; j<nrhs; j++) {
 49:     for (i=0; i<n; i++) {
 50:       PetscRandomGetValue(rand,&rval);
 51:       array[n*j+i] = rval;
 52:     }
 53:   }
 54:   MatDenseRestoreArray(RHS,&array);

 56:   MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&SOLU);

 58:   /* create matrix */
 59:   MatCreateSeqDense(PETSC_COMM_WORLD,m,n,NULL,&mat);
 60:   MatGetOwnershipRange(mat,&rstart,&rend);
 61:   for (i=rstart; i<rend; i++) {
 62:     value = (PetscReal)i+1;
 63:     MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);
 64:   }
 65:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 66:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

 68:   MatGetInfo(mat,MAT_LOCAL,&info);
 69:   PetscPrintf(PETSC_COMM_WORLD,"matrix nonzeros = %D, allocated nonzeros = %D\n",
 70:                      (PetscInt)info.nz_used,(PetscInt)info.nz_allocated);

 72:   /* Cholesky factorization - perm and factinfo are ignored by LAPACK */
 73:   /* in-place Cholesky */
 74:   MatMult(mat,x,b);
 75:   MatDuplicate(mat,MAT_COPY_VALUES,&F);
 76:   MatCholeskyFactor(F,0,0);
 77:   MatSolve(F,b,y);
 78:   MatDestroy(&F);
 79:   value = -1.0; VecAXPY(y,value,x);
 80:   VecNorm(y,NORM_2,&norm);
 81:   if (norm > tol) {
 82:     PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for Cholesky %g\n",(double)norm);
 83:   }

 85:   /* out-place Cholesky */
 86:   MatGetFactor(mat,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F);
 87:   MatCholeskyFactorSymbolic(F,mat,0,0);
 88:   MatCholeskyFactorNumeric(F,mat,0);
 89:   MatSolve(F,b,y);
 90:   value = -1.0; VecAXPY(y,value,x);
 91:   VecNorm(y,NORM_2,&norm);
 92:   if (norm > tol) {
 93:     PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for Cholesky %g\n",(double)norm);
 94:   }
 95:   MatDestroy(&F);

 97:   /* LU factorization - perms and factinfo are ignored by LAPACK */
 98:   i    = m-1; value = 1.0;
 99:   MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);
100:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
101:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
102:   MatMult(mat,x,b);
103:   MatDuplicate(mat,MAT_COPY_VALUES,&F);

105:   /* in-place LU */
106:   MatLUFactor(F,0,0,0);
107:   MatSolve(F,b,y);
108:   value = -1.0; VecAXPY(y,value,x);
109:   VecNorm(y,NORM_2,&norm);
110:   if (norm > tol) {
111:     PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for LU %g\n",(double)norm);
112:   }
113:   MatMatSolve(F,RHS,SOLU);
114:   MatDenseGetArray(SOLU,&solu_array);
115:   MatDenseGetArray(RHS,&rhs_array);
116:   for (j=0; j<nrhs; j++) {
117:     VecPlaceArray(y,solu_array+j*m);
118:     VecPlaceArray(b,rhs_array+j*m);

120:     MatMult(mat,y,ytmp);
121:     VecAXPY(ytmp,-1.0,b); /* ytmp = mat*SOLU[:,j] - RHS[:,j] */
122:     VecNorm(ytmp,NORM_2,&norm);
123:     if (norm > tol) {
124:       PetscPrintf(PETSC_COMM_WORLD,"Error: Norm of residual for LU %g\n",(double)norm);
125:     }

127:     VecResetArray(b);
128:     VecResetArray(y);
129:   }
130:   MatDenseRestoreArray(RHS,&rhs_array);
131:   MatDenseRestoreArray(SOLU,&solu_array);

133:   MatDestroy(&F);

135:   /* out-place LU */
136:   MatGetFactor(mat,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
137:   MatLUFactorSymbolic(F,mat,0,0,0);
138:   MatLUFactorNumeric(F,mat,0);
139:   MatSolve(F,b,y);
140:   value = -1.0; VecAXPY(y,value,x);
141:   VecNorm(y,NORM_2,&norm);
142:   if (norm > tol) {
143:     PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for LU %g\n",(double)norm);
144:   }

146:   /* free space */
147:   MatDestroy(&F);
148:   MatDestroy(&mat);
149:   MatDestroy(&RHS);
150:   MatDestroy(&SOLU);
151:   PetscRandomDestroy(&rand);
152:   VecDestroy(&x);
153:   VecDestroy(&b);
154:   VecDestroy(&y);
155:   VecDestroy(&ytmp);
156:   PetscFinalize();
157:   return 0;
158: }