Actual source code: ex145.c
petsc-3.10.5 2019-03-28
2: static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for an Elemental dense matrix.\n\n";
4: #include <petscmat.h>
6: int main(int argc,char **argv)
7: {
8: Mat A,F,B,X,C,Aher,G;
9: Vec b,x,c,d,e;
11: PetscInt m = 5,n,p,i,j,nrows,ncols;
12: PetscScalar *v,*barray,rval;
13: PetscReal norm,tol=1.e-11;
14: PetscMPIInt size,rank;
15: PetscRandom rand;
16: const PetscInt *rows,*cols;
17: IS isrows,iscols;
18: PetscBool mats_view=PETSC_FALSE;
19: MatFactorInfo finfo;
21: PetscInitialize(&argc,&argv,(char*) 0,help);if (ierr) return ierr;
22: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
23: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
26: PetscRandomSetFromOptions(rand);
28: /* Get local dimensions of matrices */
29: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
30: n = m;
31: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
32: p = m/2;
33: PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
34: PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view);
36: /* Create matrix A */
37: PetscPrintf(PETSC_COMM_WORLD," Create Elemental matrix A\n");
38: MatCreate(PETSC_COMM_WORLD,&A);
39: MatSetSizes(A,m,n,PETSC_DECIDE,PETSC_DECIDE);
40: MatSetType(A,MATELEMENTAL);
41: MatSetFromOptions(A);
42: MatSetUp(A);
43: /* Set local matrix entries */
44: MatGetOwnershipIS(A,&isrows,&iscols);
45: ISGetLocalSize(isrows,&nrows);
46: ISGetIndices(isrows,&rows);
47: ISGetLocalSize(iscols,&ncols);
48: ISGetIndices(iscols,&cols);
49: PetscMalloc1(nrows*ncols,&v);
50: for (i=0; i<nrows; i++) {
51: for (j=0; j<ncols; j++) {
52: PetscRandomGetValue(rand,&rval);
53: v[i*ncols+j] = rval;
54: }
55: }
56: MatSetValues(A,nrows,rows,ncols,cols,v,INSERT_VALUES);
57: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
58: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
59: ISRestoreIndices(isrows,&rows);
60: ISRestoreIndices(iscols,&cols);
61: ISDestroy(&isrows);
62: ISDestroy(&iscols);
63: PetscFree(v);
64: if (mats_view) {
65: PetscPrintf(PETSC_COMM_WORLD, "A: nrows %d, m %d; ncols %d, n %d\n",nrows,m,ncols,n);
66: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
67: }
69: /* Create rhs matrix B */
70: PetscPrintf(PETSC_COMM_WORLD," Create rhs matrix B\n");
71: MatCreate(PETSC_COMM_WORLD,&B);
72: MatSetSizes(B,m,p,PETSC_DECIDE,PETSC_DECIDE);
73: MatSetType(B,MATELEMENTAL);
74: MatSetFromOptions(B);
75: MatSetUp(B);
76: MatGetOwnershipIS(B,&isrows,&iscols);
77: ISGetLocalSize(isrows,&nrows);
78: ISGetIndices(isrows,&rows);
79: ISGetLocalSize(iscols,&ncols);
80: ISGetIndices(iscols,&cols);
81: PetscMalloc1(nrows*ncols,&v);
82: for (i=0; i<nrows; i++) {
83: for (j=0; j<ncols; j++) {
84: PetscRandomGetValue(rand,&rval);
85: v[i*ncols+j] = rval;
86: }
87: }
88: MatSetValues(B,nrows,rows,ncols,cols,v,INSERT_VALUES);
89: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
90: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
91: ISRestoreIndices(isrows,&rows);
92: ISRestoreIndices(iscols,&cols);
93: ISDestroy(&isrows);
94: ISDestroy(&iscols);
95: PetscFree(v);
96: if (mats_view) {
97: PetscPrintf(PETSC_COMM_WORLD, "B: nrows %d, m %d; ncols %d, p %d\n",nrows,m,ncols,p);
98: MatView(B,PETSC_VIEWER_STDOUT_WORLD);
99: }
101: /* Create rhs vector b and solution x (same size as b) */
102: VecCreate(PETSC_COMM_WORLD,&b);
103: VecSetSizes(b,m,PETSC_DECIDE);
104: VecSetFromOptions(b);
105: VecGetArray(b,&barray);
106: for (j=0; j<m; j++) {
107: PetscRandomGetValue(rand,&rval);
108: barray[j] = rval;
109: }
110: VecRestoreArray(b,&barray);
111: VecAssemblyBegin(b);
112: VecAssemblyEnd(b);
113: if (mats_view) {
114: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] b: m %d\n",rank,m);
115: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
116: VecView(b,PETSC_VIEWER_STDOUT_WORLD);
117: }
118: VecDuplicate(b,&x);
120: /* Create matrix X - same size as B */
121: PetscPrintf(PETSC_COMM_WORLD," Create solution matrix X\n");
122: MatCreate(PETSC_COMM_WORLD,&X);
123: MatSetSizes(X,m,p,PETSC_DECIDE,PETSC_DECIDE);
124: MatSetType(X,MATELEMENTAL);
125: MatSetFromOptions(X);
126: MatSetUp(X);
127: MatAssemblyBegin(X,MAT_FINAL_ASSEMBLY);
128: MatAssemblyEnd(X,MAT_FINAL_ASSEMBLY);
130: /* Cholesky factorization */
131: /*------------------------*/
132: PetscPrintf(PETSC_COMM_WORLD," Create Elemental matrix Aher\n");
133: MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&Aher);
134: MatAXPY(Aher,1.0,A,SAME_NONZERO_PATTERN); /* Aher = A + A^T */
135: if(!rank) { /* add 100.0 to diagonals of Aher to make it spd */
137: /* TODO: Replace this with a call to El::ShiftDiagonal( A, 100. ),
138: or at least pre-allocate the right amount of space */
139: PetscInt M,N;
140: MatGetSize(Aher,&M,&N);
141: for (i=0; i<M; i++) {
142: rval = 100.0;
143: MatSetValues(Aher,1,&i,1,&i,&rval,ADD_VALUES);
144: }
145: }
146: MatAssemblyBegin(Aher,MAT_FINAL_ASSEMBLY);
147: MatAssemblyEnd(Aher,MAT_FINAL_ASSEMBLY);
148: if (mats_view) {
149: PetscPrintf(PETSC_COMM_WORLD, "Aher:\n");
150: MatView(Aher,PETSC_VIEWER_STDOUT_WORLD);
151: }
153: /* Cholesky factorization */
154: /*------------------------*/
155: PetscPrintf(PETSC_COMM_WORLD," Test Cholesky Solver \n");
156: /* In-place Cholesky */
157: /* Create matrix factor G, then copy Aher to G */
158: MatCreate(PETSC_COMM_WORLD,&G);
159: MatSetSizes(G,m,n,PETSC_DECIDE,PETSC_DECIDE);
160: MatSetType(G,MATELEMENTAL);
161: MatSetFromOptions(G);
162: MatSetUp(G);
163: MatAssemblyBegin(G,MAT_FINAL_ASSEMBLY);
164: MatAssemblyEnd(G,MAT_FINAL_ASSEMBLY);
165: MatCopy(Aher,G,SAME_NONZERO_PATTERN);
167: /* Only G = U^T * U is implemented for now */
168: MatCholeskyFactor(G,0,0);
169: if (mats_view) {
170: PetscPrintf(PETSC_COMM_WORLD, "Cholesky Factor G:\n");
171: MatView(G,PETSC_VIEWER_STDOUT_WORLD);
172: }
174: /* Solve U^T * U x = b and U^T * U X = B */
175: MatSolve(G,b,x);
176: MatMatSolve(G,B,X);
177: MatDestroy(&G);
179: /* Out-place Cholesky */
180: MatGetFactor(Aher,MATSOLVERELEMENTAL,MAT_FACTOR_CHOLESKY,&G);
181: MatCholeskyFactorSymbolic(G,Aher,0,&finfo);
182: MatCholeskyFactorNumeric(G,Aher,&finfo);
183: if (mats_view) {
184: MatView(G,PETSC_VIEWER_STDOUT_WORLD);
185: }
186: MatSolve(G,b,x);
187: MatMatSolve(G,B,X);
188: MatDestroy(&G);
190: /* Check norm(Aher*x - b) */
191: VecCreate(PETSC_COMM_WORLD,&c);
192: VecSetSizes(c,m,PETSC_DECIDE);
193: VecSetFromOptions(c);
194: MatMult(Aher,x,c);
195: VecAXPY(c,-1.0,b);
196: VecNorm(c,NORM_1,&norm);
197: if (norm > tol) {
198: PetscPrintf(PETSC_COMM_WORLD,"Warning: |Aher*x - b| for Cholesky %g\n",(double)norm);
199: }
201: /* Check norm(Aher*X - B) */
202: MatMatMult(Aher,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
203: MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN);
204: MatNorm(C,NORM_1,&norm);
205: if (norm > tol) {
206: PetscPrintf(PETSC_COMM_WORLD,"Warning: |Aher*X - B| for Cholesky %g\n",(double)norm);
207: }
209: /* LU factorization */
210: /*------------------*/
211: PetscPrintf(PETSC_COMM_WORLD," Test LU Solver \n");
212: /* In-place LU */
213: /* Create matrix factor F, then copy A to F */
214: MatCreate(PETSC_COMM_WORLD,&F);
215: MatSetSizes(F,m,n,PETSC_DECIDE,PETSC_DECIDE);
216: MatSetType(F,MATELEMENTAL);
217: MatSetFromOptions(F);
218: MatSetUp(F);
219: MatAssemblyBegin(F,MAT_FINAL_ASSEMBLY);
220: MatAssemblyEnd(F,MAT_FINAL_ASSEMBLY);
221: MatCopy(A,F,SAME_NONZERO_PATTERN);
222: /* Create vector d to test MatSolveAdd() */
223: VecDuplicate(x,&d);
224: VecCopy(x,d);
226: /* PF=LU or F=LU factorization - perms is ignored by Elemental;
227: set finfo.dtcol !0 or 0 to enable/disable partial pivoting */
228: finfo.dtcol = 0.1;
229: MatLUFactor(F,0,0,&finfo);
231: /* Solve LUX = PB or LUX = B */
232: MatSolveAdd(F,b,d,x);
233: MatMatSolve(F,B,X);
234: MatDestroy(&F);
236: /* Check norm(A*X - B) */
237: VecCreate(PETSC_COMM_WORLD,&e);
238: VecSetSizes(e,m,PETSC_DECIDE);
239: VecSetFromOptions(e);
240: MatMult(A,x,c);
241: MatMult(A,d,e);
242: VecAXPY(c,-1.0,e);
243: VecAXPY(c,-1.0,b);
244: VecNorm(c,NORM_1,&norm);
245: if (norm > tol) {
246: PetscPrintf(PETSC_COMM_WORLD,"Warning: |A*x - b| for LU %g\n",(double)norm);
247: }
248: MatMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
249: MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN);
250: MatNorm(C,NORM_1,&norm);
251: if (norm > tol) {
252: PetscPrintf(PETSC_COMM_WORLD,"Warning: |A*X - B| for LU %g\n",(double)norm);
253: }
255: /* Out-place LU */
256: MatGetFactor(A,MATSOLVERELEMENTAL,MAT_FACTOR_LU,&F);
257: MatLUFactorSymbolic(F,A,0,0,&finfo);
258: MatLUFactorNumeric(F,A,&finfo);
259: if (mats_view) {
260: MatView(F,PETSC_VIEWER_STDOUT_WORLD);
261: }
262: MatSolve(F,b,x);
263: MatMatSolve(F,B,X);
264: MatDestroy(&F);
266: /* Free space */
267: MatDestroy(&A);
268: MatDestroy(&Aher);
269: MatDestroy(&B);
270: MatDestroy(&C);
271: MatDestroy(&X);
272: VecDestroy(&b);
273: VecDestroy(&c);
274: VecDestroy(&d);
275: VecDestroy(&e);
276: VecDestroy(&x);
277: PetscRandomDestroy(&rand);
278: PetscFinalize();
279: return ierr;
280: }
282: /*TEST
284: build:
285: requires: elemental
287: test:
288: nsize: 2
289: output_file: output/ex145.out
291: test:
292: suffix: 2
293: nsize: 6
294: output_file: output/ex145.out
296: TEST*/