Actual source code: ex39.c
2: static char help[] = "Tests Elemental interface.\n\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Mat Cdense,B,C,Ct;
9: Vec d;
10: PetscInt i,j,m = 5,n,nrows,ncols;
11: const PetscInt *rows,*cols;
12: IS isrows,iscols;
13: PetscScalar *v;
14: PetscMPIInt rank,size;
15: PetscReal Cnorm;
16: PetscBool flg,mats_view=PETSC_FALSE;
18: PetscInitialize(&argc,&args,(char*)0,help);
19: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
20: MPI_Comm_size(PETSC_COMM_WORLD,&size);
21: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
22: n = m;
23: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
24: PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view);
26: MatCreate(PETSC_COMM_WORLD,&C);
27: MatSetSizes(C,m,n,PETSC_DECIDE,PETSC_DECIDE);
28: MatSetType(C,MATELEMENTAL);
29: MatSetFromOptions(C);
30: MatSetUp(C);
31: MatGetOwnershipIS(C,&isrows,&iscols);
32: ISGetLocalSize(isrows,&nrows);
33: ISGetIndices(isrows,&rows);
34: ISGetLocalSize(iscols,&ncols);
35: ISGetIndices(iscols,&cols);
36: PetscMalloc1(nrows*ncols,&v);
37: #if defined(PETSC_USE_COMPLEX)
38: PetscRandom rand;
39: PetscScalar rval;
40: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
41: PetscRandomSetFromOptions(rand);
42: for (i=0; i<nrows; i++) {
43: for (j=0; j<ncols; j++) {
44: PetscRandomGetValue(rand,&rval);
45: v[i*ncols+j] = rval;
46: }
47: }
48: PetscRandomDestroy(&rand);
49: #else
50: for (i=0; i<nrows; i++) {
51: for (j=0; j<ncols; j++) {
52: v[i*ncols+j] = (PetscReal)(10000*rank+100*rows[i]+cols[j]);
53: }
54: }
55: #endif
56: MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES);
57: ISRestoreIndices(isrows,&rows);
58: ISRestoreIndices(iscols,&cols);
59: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
60: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
61: ISDestroy(&isrows);
62: ISDestroy(&iscols);
64: /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */
65: MatDuplicate(C,MAT_COPY_VALUES,&B);
66: if (mats_view) {
67: PetscPrintf(PETSC_COMM_WORLD,"Duplicated C:\n");
68: MatView(B,PETSC_VIEWER_STDOUT_WORLD);
69: }
70: MatDestroy(&B);
71: MatConvert(C,MATMPIDENSE,MAT_INITIAL_MATRIX,&Cdense);
72: MatMultEqual(C,Cdense,5,&flg);
75: /* Test MatNorm() */
76: MatNorm(C,NORM_1,&Cnorm);
78: /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */
79: MatTranspose(C,MAT_INITIAL_MATRIX,&Ct);
80: MatConjugate(Ct);
81: if (mats_view) {
82: PetscPrintf(PETSC_COMM_WORLD,"C's Transpose Conjugate:\n");
83: MatView(Ct,PETSC_VIEWER_STDOUT_WORLD);
84: }
85: MatZeroEntries(Ct);
86: VecCreate(PETSC_COMM_WORLD,&d);
87: VecSetSizes(d,m>n ? n : m,PETSC_DECIDE);
88: VecSetFromOptions(d);
89: MatGetDiagonal(C,d);
90: if (mats_view) {
91: PetscPrintf(PETSC_COMM_WORLD,"Diagonal of C:\n");
92: VecView(d,PETSC_VIEWER_STDOUT_WORLD);
93: }
94: if (m>n) {
95: MatDiagonalScale(C,NULL,d);
96: } else {
97: MatDiagonalScale(C,d,NULL);
98: }
99: if (mats_view) {
100: PetscPrintf(PETSC_COMM_WORLD,"Diagonal Scaled C:\n");
101: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
102: }
104: /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
105: MatCreate(PETSC_COMM_WORLD,&B);
106: MatSetSizes(B,m,n,PETSC_DECIDE,PETSC_DECIDE);
107: MatSetType(B,MATELEMENTAL);
108: MatSetFromOptions(B);
109: MatSetUp(B);
110: MatGetOwnershipIS(B,&isrows,&iscols);
111: ISGetLocalSize(isrows,&nrows);
112: ISGetIndices(isrows,&rows);
113: ISGetLocalSize(iscols,&ncols);
114: ISGetIndices(iscols,&cols);
115: for (i=0; i<nrows; i++) {
116: for (j=0; j<ncols; j++) {
117: v[i*ncols+j] = (PetscReal)(1000*rows[i]+cols[j]);
118: }
119: }
120: MatSetValues(B,nrows,rows,ncols,cols,v,INSERT_VALUES);
121: PetscFree(v);
122: ISRestoreIndices(isrows,&rows);
123: ISRestoreIndices(iscols,&cols);
124: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
125: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
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);
144: if (mats_view) {
145: PetscPrintf(PETSC_COMM_WORLD,"C MatMatTransposeMult C:\n");
146: MatView(B,PETSC_VIEWER_STDOUT_WORLD);
147: }
149: MatDestroy(&Cdense);
150: PetscFree(v);
151: MatDestroy(&B);
152: MatDestroy(&C);
153: MatDestroy(&Ct);
154: VecDestroy(&d);
155: PetscFinalize();
156: return 0;
157: }
159: /*TEST
161: test:
162: nsize: 2
163: args: -m 3 -n 2
164: requires: elemental
166: test:
167: suffix: 2
168: nsize: 6
169: args: -m 2 -n 3
170: requires: elemental
172: test:
173: suffix: 3
174: nsize: 1
175: args: -m 2 -n 3
176: requires: elemental
177: output_file: output/ex39_1.out
179: TEST*/