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>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

131:   MatDestroy(&F);

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

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