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