Actual source code: ex242.c
2: static char help[] = "Tests ScaLAPACK interface.\n\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Mat Cdense,Caij,B,C,Ct,Asub;
9: Vec d;
10: PetscInt i,j,M = 5,N,mb = 2,nb,nrows,ncols,mloc,nloc;
11: const PetscInt *rows,*cols;
12: IS isrows,iscols;
13: PetscScalar *v;
14: PetscMPIInt rank,color;
15: PetscReal Cnorm;
16: PetscBool flg,mats_view=PETSC_FALSE;
17: MPI_Comm subcomm;
19: PetscInitialize(&argc,&args,(char*)0,help);
20: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
21: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
22: N = M;
23: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
24: PetscOptionsGetInt(NULL,NULL,"-mb",&mb,NULL);
25: nb = mb;
26: PetscOptionsGetInt(NULL,NULL,"-nb",&nb,NULL);
27: PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view);
29: MatCreate(PETSC_COMM_WORLD,&C);
30: MatSetType(C,MATSCALAPACK);
31: mloc = PETSC_DECIDE;
32: PetscSplitOwnershipEqual(PETSC_COMM_WORLD,&mloc,&M);
33: nloc = PETSC_DECIDE;
34: PetscSplitOwnershipEqual(PETSC_COMM_WORLD,&nloc,&N);
35: MatSetSizes(C,mloc,nloc,M,N);
36: MatScaLAPACKSetBlockSizes(C,mb,nb);
37: MatSetFromOptions(C);
38: MatSetUp(C);
39: /*MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&C); */
41: MatGetOwnershipIS(C,&isrows,&iscols);
42: ISGetLocalSize(isrows,&nrows);
43: ISGetIndices(isrows,&rows);
44: ISGetLocalSize(iscols,&ncols);
45: ISGetIndices(iscols,&cols);
46: PetscMalloc1(nrows*ncols,&v);
47: for (i=0;i<nrows;i++) {
48: for (j=0;j<ncols;j++) v[i*ncols+j] = (PetscReal)(rows[i]+1+(cols[j]+1)*0.01);
49: }
50: MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES);
51: PetscFree(v);
52: ISRestoreIndices(isrows,&rows);
53: ISRestoreIndices(iscols,&cols);
54: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
55: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
56: ISDestroy(&isrows);
57: ISDestroy(&iscols);
59: /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */
60: MatDuplicate(C,MAT_COPY_VALUES,&B);
61: if (mats_view) {
62: PetscPrintf(PETSC_COMM_WORLD,"Duplicated C:\n");
63: MatView(B,PETSC_VIEWER_STDOUT_WORLD);
64: }
65: MatDestroy(&B);
66: MatConvert(C,MATDENSE,MAT_INITIAL_MATRIX,&Cdense);
67: MatMultEqual(C,Cdense,5,&flg);
70: /* Test MatNorm() */
71: MatNorm(C,NORM_1,&Cnorm);
73: /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */
74: MatTranspose(C,MAT_INITIAL_MATRIX,&Ct);
75: MatConjugate(Ct);
76: if (mats_view) {
77: PetscPrintf(PETSC_COMM_WORLD,"C's Transpose Conjugate:\n");
78: MatView(Ct,PETSC_VIEWER_STDOUT_WORLD);
79: }
80: MatZeroEntries(Ct);
81: if (M>N) MatCreateVecs(C,&d,NULL);
82: else MatCreateVecs(C,NULL,&d);
83: MatGetDiagonal(C,d);
84: if (mats_view) {
85: PetscPrintf(PETSC_COMM_WORLD,"Diagonal of C:\n");
86: VecView(d,PETSC_VIEWER_STDOUT_WORLD);
87: }
88: if (M>N) {
89: MatDiagonalScale(C,NULL,d);
90: } else {
91: MatDiagonalScale(C,d,NULL);
92: }
93: if (mats_view) {
94: PetscPrintf(PETSC_COMM_WORLD,"Diagonal Scaled C:\n");
95: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
96: }
98: /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
99: MatCreate(PETSC_COMM_WORLD,&B);
100: MatSetType(B,MATSCALAPACK);
101: MatSetSizes(B,mloc,nloc,PETSC_DECIDE,PETSC_DECIDE);
102: MatScaLAPACKSetBlockSizes(B,mb,nb);
103: MatSetFromOptions(B);
104: MatSetUp(B);
105: /* MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&B); */
106: MatGetOwnershipIS(B,&isrows,&iscols);
107: ISGetLocalSize(isrows,&nrows);
108: ISGetIndices(isrows,&rows);
109: ISGetLocalSize(iscols,&ncols);
110: ISGetIndices(iscols,&cols);
111: PetscMalloc1(nrows*ncols,&v);
112: for (i=0;i<nrows;i++) {
113: for (j=0;j<ncols;j++) v[i*ncols+j] = (PetscReal)(1000*rows[i]+cols[j]);
114: }
115: MatSetValues(B,nrows,rows,ncols,cols,v,INSERT_VALUES);
116: PetscFree(v);
117: ISRestoreIndices(isrows,&rows);
118: ISRestoreIndices(iscols,&cols);
119: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
120: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
121: if (mats_view) {
122: PetscPrintf(PETSC_COMM_WORLD,"B:\n");
123: MatView(B,PETSC_VIEWER_STDOUT_WORLD);
124: }
125: MatAXPY(B,2.5,C,SAME_NONZERO_PATTERN);
126: MatAYPX(B,3.75,C,SAME_NONZERO_PATTERN);
127: MatConvert(B,MATDENSE,MAT_INPLACE_MATRIX,&B);
128: if (mats_view) {
129: PetscPrintf(PETSC_COMM_WORLD,"B after MatAXPY and MatAYPX:\n");
130: MatView(B,PETSC_VIEWER_STDOUT_WORLD);
131: }
132: ISDestroy(&isrows);
133: ISDestroy(&iscols);
134: MatDestroy(&B);
136: /* Test MatMatTransposeMult(): B = C*C^T */
137: MatMatTransposeMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B);
138: MatScale(C,2.0);
139: MatMatTransposeMult(C,C,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);
140: MatMatTransposeMultEqual(C,C,B,10,&flg);
143: if (mats_view) {
144: PetscPrintf(PETSC_COMM_WORLD,"C MatMatTransposeMult C:\n");
145: MatView(B,PETSC_VIEWER_STDOUT_WORLD);
146: }
148: /* Test MatMult() */
149: MatComputeOperator(C,MATAIJ,&Caij);
150: MatMultEqual(C,Caij,5,&flg);
152: MatMultTransposeEqual(C,Caij,5,&flg);
155: /* Test MatMultAdd() and MatMultTransposeAddEqual() */
156: MatMultAddEqual(C,Caij,5,&flg);
158: MatMultTransposeAddEqual(C,Caij,5,&flg);
161: /* Test MatMatMult() */
162: PetscOptionsHasName(NULL,NULL,"-test_matmatmult",&flg);
163: if (flg) {
164: Mat CC,CCaij;
165: MatMatMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CC);
166: MatMatMult(Caij,Caij,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCaij);
167: MatMultEqual(CC,CCaij,5,&flg);
169: MatDestroy(&CCaij);
170: MatDestroy(&CC);
171: }
173: /* Test MatCreate() on subcomm */
174: color = rank%2;
175: MPI_Comm_split(PETSC_COMM_WORLD,color,0,&subcomm);
176: if (color==0) {
177: MatCreate(subcomm,&Asub);
178: MatSetType(Asub,MATSCALAPACK);
179: mloc = PETSC_DECIDE;
180: PetscSplitOwnershipEqual(subcomm,&mloc,&M);
181: nloc = PETSC_DECIDE;
182: PetscSplitOwnershipEqual(subcomm,&nloc,&N);
183: MatSetSizes(Asub,mloc,nloc,M,N);
184: MatScaLAPACKSetBlockSizes(Asub,mb,nb);
185: MatSetFromOptions(Asub);
186: MatSetUp(Asub);
187: MatDestroy(&Asub);
188: }
190: MatDestroy(&Cdense);
191: MatDestroy(&Caij);
192: MatDestroy(&B);
193: MatDestroy(&C);
194: MatDestroy(&Ct);
195: VecDestroy(&d);
196: MPI_Comm_free(&subcomm);
197: PetscFinalize();
198: return 0;
199: }
201: /*TEST
203: build:
204: requires: scalapack
206: test:
207: nsize: 2
208: args: -mb 5 -nb 5 -M 12 -N 10
209: requires: scalapack
211: test:
212: suffix: 2
213: nsize: 6
214: args: -mb 8 -nb 6 -M 20 -N 50
215: requires: scalapack
216: output_file: output/ex242_1.out
218: test:
219: suffix: 3
220: nsize: 3
221: args: -mb 2 -nb 2 -M 20 -N 20 -test_matmatmult
222: requires: scalapack
223: output_file: output/ex242_1.out
225: TEST*/