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:   factinfo.fill = 1.0;
 48:   MatMult(mat,x,b);
 49:   MatConvert(mat,MATSAME,MAT_INITIAL_MATRIX,&fact);
 50:   MatCholeskyFactor(fact,perm,&factinfo);
 51:   MatSolve(fact,b,y);
 52:   MatDestroy(fact);
 53:   value = -1.0; VecAXPY(y,value,x);
 54:   VecNorm(y,NORM_2,&norm);
 55:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error for Cholesky %A\n",norm);
 56:   MatCholeskyFactorSymbolic(mat,perm,&factinfo,&fact);
 57:   MatCholeskyFactorNumeric(mat,&factinfo,&fact);
 58:   MatSolve(fact,b,y);
 59:   value = -1.0; VecAXPY(y,value,x);
 60:   VecNorm(y,NORM_2,&norm);
 61:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error for Cholesky %A\n",norm);
 62:   MatDestroy(fact);

 64:   i = m-1; value = 1.0;
 65:   MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);
 66:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 67:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
 68:   MatMult(mat,x,b);
 69:   MatConvert(mat,MATSAME,MAT_INITIAL_MATRIX,&fact);
 70:   luinfo.fill           = 5.0;
 71:   luinfo.dtcol          = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
 72:   luinfo.shiftnz        = 0.0;
 73:   luinfo.zeropivot      = 1.e-12;
 74:   luinfo.pivotinblocks  = 1.0;
 75:   luinfo.shiftpd        = 0.0; /* false */
 76:   luinfo.shift_fraction = 0.0;
 77:   MatLUFactor(fact,perm,perm,&luinfo);
 78:   MatSolve(fact,b,y);
 79:   value = -1.0; VecAXPY(y,value,x);
 80:   VecNorm(y,NORM_2,&norm);
 81:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
 82:   MatDestroy(fact);

 84:   luinfo.fill      = 2.0;
 85:   luinfo.dtcol     = 0.0;
 86:   luinfo.shiftnz   = 0.0;
 87:   luinfo.zeropivot = 1.e-14;
 88:   luinfo.pivotinblocks = 1.0;
 89:   MatLUFactorSymbolic(mat,perm,perm,&luinfo,&fact);
 90:   MatLUFactorNumeric(mat,&luinfo,&fact);
 91:   MatSolve(fact,b,y);
 92:   value = -1.0; VecAXPY(y,value,x);
 93:   VecNorm(y,NORM_2,&norm);
 94:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
 95:   MatDestroy(fact);
 96:   MatDestroy(mat);
 97:   VecDestroy(x);
 98:   VecDestroy(b);
 99:   VecDestroy(y);
100:   ISDestroy(perm);

102:   PetscFinalize();
103:   return 0;
104: }
105: