1*c4762a1bSJed Brown 2*c4762a1bSJed Brown static char help[] = "Tests Cholesky factorization for a SBAIJ matrix, (bs=2).\n"; 3*c4762a1bSJed Brown /* 4*c4762a1bSJed Brown This code is modified from the code contributed by JUNWANG@uwm.edu on Apr 13, 2007 5*c4762a1bSJed Brown */ 6*c4762a1bSJed Brown 7*c4762a1bSJed Brown #include <petscmat.h> 8*c4762a1bSJed Brown 9*c4762a1bSJed Brown int main(int argc,char **args) 10*c4762a1bSJed Brown { 11*c4762a1bSJed Brown PetscErrorCode ierr; 12*c4762a1bSJed Brown Mat mat,fact,B; 13*c4762a1bSJed Brown PetscInt ind1[2],ind2[2]; 14*c4762a1bSJed Brown PetscScalar temp[4]; 15*c4762a1bSJed Brown PetscInt nnz[3]; 16*c4762a1bSJed Brown IS perm,colp; 17*c4762a1bSJed Brown MatFactorInfo info; 18*c4762a1bSJed Brown PetscMPIInt size; 19*c4762a1bSJed Brown 20*c4762a1bSJed Brown ierr = PetscInitialize(&argc,&args,0,help);if (ierr) return ierr; 21*c4762a1bSJed Brown ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 22*c4762a1bSJed Brown if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!"); 23*c4762a1bSJed Brown 24*c4762a1bSJed Brown nnz[0]=2;nnz[1]=1;nnz[2]=1; 25*c4762a1bSJed Brown ierr = MatCreateSeqSBAIJ(PETSC_COMM_SELF,2,6,6,0,nnz,&mat);CHKERRQ(ierr); 26*c4762a1bSJed Brown 27*c4762a1bSJed Brown ind1[0]=0;ind1[1]=1; 28*c4762a1bSJed Brown temp[0]=3;temp[1]=2;temp[2]=0;temp[3]=3; 29*c4762a1bSJed Brown ierr = MatSetValues(mat,2,ind1,2,ind1,temp,INSERT_VALUES);CHKERRQ(ierr); 30*c4762a1bSJed Brown ind2[0]=4;ind2[1]=5; 31*c4762a1bSJed Brown temp[0]=1;temp[1]=1;temp[2]=2;temp[3]=1; 32*c4762a1bSJed Brown ierr = MatSetValues(mat,2,ind1,2,ind2,temp,INSERT_VALUES);CHKERRQ(ierr); 33*c4762a1bSJed Brown ind1[0]=2;ind1[1]=3; 34*c4762a1bSJed Brown temp[0]=4;temp[1]=1;temp[2]=1;temp[3]=5; 35*c4762a1bSJed Brown ierr = MatSetValues(mat,2,ind1,2,ind1,temp,INSERT_VALUES);CHKERRQ(ierr); 36*c4762a1bSJed Brown ind1[0]=4;ind1[1]=5; 37*c4762a1bSJed Brown temp[0]=5;temp[1]=1;temp[2]=1;temp[3]=6; 38*c4762a1bSJed Brown ierr = MatSetValues(mat,2,ind1,2,ind1,temp,INSERT_VALUES);CHKERRQ(ierr); 39*c4762a1bSJed Brown 40*c4762a1bSJed Brown ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 41*c4762a1bSJed Brown ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 42*c4762a1bSJed Brown 43*c4762a1bSJed Brown ierr = MatDuplicate(mat,MAT_SHARE_NONZERO_PATTERN,&B);CHKERRQ(ierr); 44*c4762a1bSJed Brown ind1[0]=0;ind1[1]=1; 45*c4762a1bSJed Brown temp[0]=3;temp[1]=2;temp[2]=0;temp[3]=3; 46*c4762a1bSJed Brown ierr = MatSetValues(mat,2,ind1,2,ind1,temp,INSERT_VALUES);CHKERRQ(ierr); 47*c4762a1bSJed Brown ind2[0]=4;ind2[1]=5; 48*c4762a1bSJed Brown temp[0]=1;temp[1]=1;temp[2]=2;temp[3]=1; 49*c4762a1bSJed Brown ierr = MatSetValues(mat,2,ind1,2,ind2,temp,INSERT_VALUES);CHKERRQ(ierr); 50*c4762a1bSJed Brown ind1[0]=2;ind1[1]=3; 51*c4762a1bSJed Brown temp[0]=4;temp[1]=1;temp[2]=1;temp[3]=5; 52*c4762a1bSJed Brown ierr = MatSetValues(mat,2,ind1,2,ind1,temp,INSERT_VALUES);CHKERRQ(ierr); 53*c4762a1bSJed Brown ind1[0]=4;ind1[1]=5; 54*c4762a1bSJed Brown temp[0]=5;temp[1]=1;temp[2]=1;temp[3]=6; 55*c4762a1bSJed Brown ierr = MatSetValues(mat,2,ind1,2,ind1,temp,INSERT_VALUES);CHKERRQ(ierr); 56*c4762a1bSJed Brown 57*c4762a1bSJed Brown ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 58*c4762a1bSJed Brown ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 59*c4762a1bSJed Brown 60*c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_WORLD,"mat: \n");CHKERRQ(ierr); 61*c4762a1bSJed Brown ierr = MatView(mat,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 62*c4762a1bSJed Brown 63*c4762a1bSJed Brown /* begin cholesky factorization */ 64*c4762a1bSJed Brown ierr = MatGetOrdering(mat,MATORDERINGNATURAL,&perm,&colp);CHKERRQ(ierr); 65*c4762a1bSJed Brown ierr = ISDestroy(&colp);CHKERRQ(ierr); 66*c4762a1bSJed Brown 67*c4762a1bSJed Brown info.fill=1.0; 68*c4762a1bSJed Brown ierr = MatGetFactor(mat,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&fact);CHKERRQ(ierr); 69*c4762a1bSJed Brown ierr = MatCholeskyFactorSymbolic(fact,mat,perm,&info);CHKERRQ(ierr); 70*c4762a1bSJed Brown ierr = MatCholeskyFactorNumeric(fact,mat,&info);CHKERRQ(ierr); 71*c4762a1bSJed Brown 72*c4762a1bSJed Brown ierr = ISDestroy(&perm);CHKERRQ(ierr); 73*c4762a1bSJed Brown ierr = MatDestroy(&mat);CHKERRQ(ierr); 74*c4762a1bSJed Brown ierr = MatDestroy(&fact);CHKERRQ(ierr); 75*c4762a1bSJed Brown ierr = MatDestroy(&B);CHKERRQ(ierr); 76*c4762a1bSJed Brown ierr = PetscFinalize(); 77*c4762a1bSJed Brown return ierr; 78*c4762a1bSJed Brown } 79