Actual source code: ex1.c

  2: static char help[] = "Tests LU and Cholesky factorization for a dense matrix.\n\n";

 4:  #include petscmat.h

  8: int main(int argc,char **argv)
  9: {
 10:   Mat            mat,fact;
 11:   MatInfo        info;
 13:   PetscInt       m = 10,n = 10,i = 4,rstart,rend;
 14:   PetscScalar    value = 1.0;
 15:   Vec            x,y,b;
 16:   PetscReal      norm;
 17:   IS             perm;
 18:   MatFactorInfo  luinfo,factinfo;

 20:   PetscInitialize(&argc,&argv,(char*) 0,help);

 22:   VecCreate(PETSC_COMM_WORLD,&y);
 23:   VecSetSizes(y,PETSC_DECIDE,m);
 24:   VecSetFromOptions(y);
 25:   VecDuplicate(y,&x);
 26:   VecSet(x,value);
 27:   VecCreate(PETSC_COMM_WORLD,&b);
 28:   VecSetSizes(b,PETSC_DECIDE,n);
 29:   VecSetFromOptions(b);
 30:   ISCreateStride(PETSC_COMM_WORLD,m,0,1,&perm);

 32:   MatCreateSeqDense(PETSC_COMM_WORLD,m,n,PETSC_NULL,&mat);

 34:   MatGetOwnershipRange(mat,&rstart,&rend);
 35:   for (i=rstart; i<rend; i++) {
 36:     value = (PetscReal)i+1;
 37:     MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);
 38:   }
 39:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 40:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

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

 46:   /* Cholesky factorization is not yet in place for this matrix format */
 47:   MatFactorInfoInitialize(&factinfo);
 48:   factinfo.fill = 1.0;
 49:   MatMult(mat,x,b);
 50:   MatConvert(mat,MATSAME,MAT_INITIAL_MATRIX,&fact);
 51:   MatCholeskyFactor(fact,perm,&factinfo);
 52:   MatSolve(fact,b,y);
 53:   MatDestroy(fact);
 54:   value = -1.0; VecAXPY(y,value,x);
 55:   VecNorm(y,NORM_2,&norm);
 56:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error for Cholesky %A\n",norm);
 57:   MatGetFactor(mat,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&fact);
 58:   MatCholeskyFactorSymbolic(fact,mat,perm,&factinfo);
 59:   MatCholeskyFactorNumeric(fact,mat,&factinfo);
 60:   MatSolve(fact,b,y);
 61:   value = -1.0; VecAXPY(y,value,x);
 62:   VecNorm(y,NORM_2,&norm);
 63:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error for Cholesky %A\n",norm);
 64:   MatDestroy(fact);

 66:   i = m-1; value = 1.0;
 67:   MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);
 68:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 69:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
 70:   MatMult(mat,x,b);
 71:   MatConvert(mat,MATSAME,MAT_INITIAL_MATRIX,&fact);

 73:   MatFactorInfoInitialize(&luinfo);
 74:   luinfo.fill           = 1.0;
 75:   luinfo.dtcol          = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
 76:   luinfo.zeropivot      = 1.e-12;
 77:   luinfo.shifttype      = (PetscReal)MAT_SHIFT_INBLOCKS;
 78:   luinfo.shiftamount    = 1.e-12;
 79: 
 80:   MatLUFactor(fact,perm,perm,&luinfo);
 81:   MatSolve(fact,b,y);
 82:   value = -1.0; VecAXPY(y,value,x);
 83:   VecNorm(y,NORM_2,&norm);
 84:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
 85:   MatDestroy(fact);

 87:   luinfo.fill        = 1.0;
 88:   luinfo.dtcol       = 0.0;
 89:   luinfo.zeropivot   = 1.e-14;
 90:   luinfo.shifttype   = (PetscReal)MAT_SHIFT_INBLOCKS;
 91:   luinfo.shiftamount = 1.e-12;
 92:   MatGetFactor(mat,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&fact);
 93:   MatLUFactorSymbolic(fact,mat,perm,perm,&luinfo);
 94:   MatLUFactorNumeric(fact,mat,&luinfo);
 95:   MatSolve(fact,b,y);
 96:   value = -1.0; VecAXPY(y,value,x);
 97:   VecNorm(y,NORM_2,&norm);
 98:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
 99:   MatDestroy(fact);
100:   MatDestroy(mat);
101:   VecDestroy(x);
102:   VecDestroy(b);
103:   VecDestroy(y);
104:   ISDestroy(perm);
105:   PetscFinalize();
106:   return 0;
107: }
108: