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