Actual source code: ex127.c
1: static char help[] = "Test MatMult() for Hermitian matrix.\n\n";
3: #include <petscmat.h>
5: int main(int argc,char **args)
6: {
7: Mat A,As;
8: PetscBool flg;
9: PetscMPIInt size;
10: PetscInt i,j;
11: PetscScalar v,sigma2;
12: PetscReal h2,sigma1=100.0;
13: PetscInt dim,Ii,J,n = 3,rstart,rend;
15: PetscInitialize(&argc,&args,(char*)0,help);
16: MPI_Comm_size(PETSC_COMM_WORLD,&size);
17: PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);
18: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
19: dim = n*n;
21: MatCreate(PETSC_COMM_WORLD,&A);
22: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
23: MatSetType(A,MATAIJ);
24: MatSetFromOptions(A);
25: MatSetUp(A);
27: sigma2 = 10.0*PETSC_i;
28: h2 = 1.0/((n+1)*(n+1));
30: MatGetOwnershipRange(A,&rstart,&rend);
31: for (Ii=rstart; Ii<rend; Ii++) {
32: v = -1.0; i = Ii/n; j = Ii - i*n;
33: if (i>0) {
34: J = Ii-n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
35: }
36: if (i<n-1) {
37: J = Ii+n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
38: }
39: if (j>0) {
40: J = Ii-1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
41: }
42: if (j<n-1) {
43: J = Ii+1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
44: }
45: v = 4.0 - sigma1*h2;
46: MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
47: }
48: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
49: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
51: /* Check whether A is symmetric */
52: PetscOptionsHasName(NULL,NULL, "-check_symmetric", &flg);
53: if (flg) {
54: MatIsSymmetric(A,0.0,&flg);
56: }
57: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
59: /* make A complex Hermitian */
60: Ii = 0; J = dim-1;
61: if (Ii >= rstart && Ii < rend) {
62: v = sigma2*h2; /* RealPart(v) = 0.0 */
63: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
64: v = -sigma2*h2;
65: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
66: }
68: Ii = dim-2; J = dim-1;
69: if (Ii >= rstart && Ii < rend) {
70: v = sigma2*h2; /* RealPart(v) = 0.0 */
71: MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);
72: v = -sigma2*h2;
73: MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);
74: }
76: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
77: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
78: MatViewFromOptions(A,NULL,"-disp_mat");
80: /* Check whether A is Hermitian, then set A->hermitian flag */
81: PetscOptionsHasName(NULL,NULL, "-check_Hermitian", &flg);
82: if (flg && size == 1) {
83: MatIsHermitian(A,0.0,&flg);
85: }
86: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
88: #if defined(PETSC_HAVE_SUPERLU_DIST)
89: /* Test Cholesky factorization */
90: PetscOptionsHasName(NULL,NULL, "-test_choleskyfactor", &flg);
91: if (flg) {
92: Mat F;
93: IS perm,iperm;
94: MatFactorInfo info;
95: PetscInt nneg,nzero,npos;
97: MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_CHOLESKY,&F);
98: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
99: MatCholeskyFactorSymbolic(F,A,perm,&info);
100: MatCholeskyFactorNumeric(F,A,&info);
102: MatGetInertia(F,&nneg,&nzero,&npos);
103: PetscPrintf(PETSC_COMM_WORLD," MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n",nneg,nzero,npos);
104: MatDestroy(&F);
105: ISDestroy(&perm);
106: ISDestroy(&iperm);
107: }
108: #endif
110: /* Create a Hermitian matrix As in sbaij format */
111: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&As);
112: MatViewFromOptions(As,NULL,"-disp_mat");
114: /* Test MatMult */
115: MatMultEqual(A,As,10,&flg);
117: MatMultAddEqual(A,As,10,&flg);
120: /* Free spaces */
121: MatDestroy(&A);
122: MatDestroy(&As);
123: PetscFinalize();
124: return 0;
125: }
127: /*TEST
129: build:
130: requires: complex
132: test:
133: args: -n 1000
134: output_file: output/ex127.out
136: test:
137: suffix: 2
138: nsize: 3
139: args: -n 1000
140: output_file: output/ex127.out
142: test:
143: suffix: superlu_dist
144: nsize: 3
145: requires: superlu_dist
146: args: -test_choleskyfactor -mat_superlu_dist_rowperm NOROWPERM
147: output_file: output/ex127_superlu_dist.out
148: TEST*/