Actual source code: ex1.c

  1: /*$Id: ex1.c,v 1.23 2001/08/07 21:30:08 bsmith Exp $*/

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

 5:  #include petscmat.h

  9: int main(int argc,char **argv)
 10: {
 11:   Mat          mat,fact;
 12:   MatInfo      info;
 13:   int          m = 10,n = 10,i = 4,ierr,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(&value,x);
 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:     (int)info.nz_used,(int)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,&fact);
 50:   MatCholeskyFactor(fact,perm,&factinfo);
 51:   MatSolve(fact,b,y);
 52:   MatDestroy(fact);
 53:   value = -1.0; VecAXPY(&value,x,y);
 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,&fact);
 58:   MatSolve(fact,b,y);
 59:   value = -1.0; VecAXPY(&value,x,y);
 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,&fact);
 70:   MatLUFactor(fact,perm,perm,PETSC_NULL);
 71:   MatSolve(fact,b,y);
 72:   value = -1.0; VecAXPY(&value,x,y);
 73:   VecNorm(y,NORM_2,&norm);
 74:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
 75:   MatDestroy(fact);

 77:   luinfo.fill = 2.0;
 78:   luinfo.dtcol = 0.0;
 79:   luinfo.damping = 0.0;
 80:   luinfo.zeropivot = 1.e-14;
 81:   luinfo.pivotinblocks = 1.0;
 82:   MatLUFactorSymbolic(mat,perm,perm,&luinfo,&fact);
 83:   MatLUFactorNumeric(mat,&fact);
 84:   MatSolve(fact,b,y);
 85:   value = -1.0; VecAXPY(&value,x,y);
 86:   VecNorm(y,NORM_2,&norm);
 87:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
 88:   MatDestroy(fact);
 89:   MatDestroy(mat);
 90:   VecDestroy(x);
 91:   VecDestroy(b);
 92:   VecDestroy(y);
 93:   ISDestroy(perm);

 95:   PetscFinalize();
 96:   return 0;
 97: }
 98: