Actual source code: ex174.cxx
petsc-3.6.1 2015-08-06
2: static char help[] = "Tests MatConvert(), MatLoad(), MatElementalHermitianGenDefEig() for MATELEMENTAL interface.\n\n";
3: /*
4: Example:
5: mpiexec -n <np> ./ex173 -fA <A_data> -fB <B_data> -vl <vl> -vu <vu> -orig_mat_type <type> -orig_mat_type <mat_type>
6: */
7:
8: #include <petscmat.h>
9: #include <petscmatelemental.h>
13: int main(int argc,char **args)
14: {
15: Mat A,Ae,B,Be,X,Xe,We,C1,C2,EVAL;
16: Vec eval;
18: PetscViewer view;
19: char file[2][PETSC_MAX_PATH_LEN];
20: PetscBool flg,flgB,isElemental,isDense,isAij,isSbaij;
21: PetscScalar one = 1.0,*Earray;
22: PetscMPIInt rank,size;
23: PetscReal vl,vu,norm;
24: PetscInt M,N,m;
26: /* Below are Elemental data types, see <elemental.hpp> */
27: El::Pencil eigtype = El::AXBX;
28: El::UpperOrLower uplo = El::UPPER;
29: El::SortType sort = El::UNSORTED; /* UNSORTED, DESCENDING, ASCENDING */
30: El::HermitianEigSubset<PetscElemScalar> subset;
31: const El::HermitianEigCtrl<PetscElemScalar> ctrl;
33: PetscInitialize(&argc,&args,(char*)0,help);
34: #if !defined(PETSC_HAVE_ELEMENTAL)
35: SETERRQ(PETSC_COMM_WORLD,1,"This example requires ELEMENTAL");
36: #endif
37: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
38: MPI_Comm_size(PETSC_COMM_WORLD,&size);
40: /* Load PETSc matrices */
41: PetscOptionsGetString(NULL,"-fA",file[0],PETSC_MAX_PATH_LEN,NULL);
42: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&view);
43: MatCreate(PETSC_COMM_WORLD,&A);
44: MatSetOptionsPrefix(A,"orig_");
45: MatSetType(A,MATAIJ);
46: MatSetFromOptions(A);
47: MatLoad(A,view);
48: PetscViewerDestroy(&view);
50: PetscOptionsGetString(NULL,"-fB",file[1],PETSC_MAX_PATH_LEN,&flgB);
51: if (flgB) {
52: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&view);
53: MatCreate(PETSC_COMM_WORLD,&B);
54: MatSetOptionsPrefix(B,"orig_");
55: MatSetType(B,MATAIJ);
56: MatSetFromOptions(B);
57: MatLoad(B,view);
58: PetscViewerDestroy(&view);
59: } else {
60: /* Create matrix B = I */
61: PetscInt rstart,rend,i;
62: MatGetSize(A,&M,&N);
63: MatGetOwnershipRange(A,&rstart,&rend);
65: MatCreate(PETSC_COMM_WORLD,&B);
66: MatSetOptionsPrefix(B,"orig_");
67: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);
68: MatSetType(B,MATAIJ);
69: MatSetFromOptions(B);
70: MatSetUp(B);
71: for (i=rstart; i<rend; i++) {
72: MatSetValues(B,1,&i,1,&i,&one,ADD_VALUES);
73: }
74: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
75: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
76: }
78: PetscObjectTypeCompare((PetscObject)A,MATELEMENTAL,&isElemental);
79: if (isElemental) {
80: Ae = A;
81: Be = B;
82: isDense = isAij = isSbaij = PETSC_FALSE;
83: } else { /* Convert AIJ/DENSE/SBAIJ matrices into Elemental matrices */
84: if (size == 1) {
85: PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&isDense);
86: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isAij);
87: PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSbaij);
88: } else {
89: PetscObjectTypeCompare((PetscObject)A,MATMPIDENSE,&isDense);
90: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isAij);
91: PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&isSbaij);
92: }
94: if (!rank) {
95: if (isDense) {
96: printf(" Convert DENSE matrices A and B into Elemental matrix... \n");
97: } else if (isAij) {
98: printf(" Convert AIJ matrices A and B into Elemental matrix... \n");
99: } else if (isSbaij) {
100: printf(" Convert SBAIJ matrices A and B into Elemental matrix... \n");
101: } else {
102: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not supported yet");
103: }
104: }
105: MatConvert(A, MATELEMENTAL, MAT_INITIAL_MATRIX, &Ae);
106: MatConvert(B, MATELEMENTAL, MAT_INITIAL_MATRIX, &Be);
108: /* Test accuracy */
109: MatMultEqual(A,Ae,5,&flg);
110: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"A != A_elemental.");
111: MatMultEqual(B,Be,5,&flg);
112: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"B != B_elemental.");
113: }
115: /* Test MatElementalHermitianGenDefEig() */
116: if (!rank) printf(" Compute Ax = lambda Bx... \n");
117: vl = -0.8, vu = -0.7;
118: PetscOptionsGetReal(NULL,"-vl",&vl,NULL);
119: PetscOptionsGetReal(NULL,"-vu",&vu,NULL);
120: subset.rangeSubset = PETSC_TRUE;
121: subset.lowerBound = vl;
122: subset.upperBound = vu;
123: MatElementalHermitianGenDefEig(eigtype,uplo,Ae,Be,&We,&Xe,sort,subset,ctrl);
124: //MatView(We,PETSC_VIEWER_STDOUT_WORLD);
126: /* Check || A*X - B*X*We || */
127: if (isAij) {
128: if (!rank) printf(" Convert Elemental matrices We and Xe into MATDENSE matrices... \n");
129: MatConvert(We,MATDENSE,MAT_INITIAL_MATRIX,&EVAL); /* EVAL is a Mx1 matrix */
130: MatConvert(Xe,MATDENSE,MAT_INITIAL_MATRIX,&X);
131: //MatView(EVAL,PETSC_VIEWER_STDOUT_WORLD);
133: MatMatMult(A,X,MAT_INITIAL_MATRIX,1.0,&C1); /* C1 = A*X */
134: MatMatMult(B,X,MAT_INITIAL_MATRIX,1.0,&C2); /* C2 = B*X */
135:
136: /* Get vector eval from matrix EVAL for MatDiagonalScale() */
137: MatGetLocalSize(EVAL,&m,NULL);
138: MatDenseGetArray(EVAL,&Earray);
139: VecCreateMPIWithArray(PetscObjectComm((PetscObject)EVAL),1,m,PETSC_DECIDE,Earray,&eval);
140: MatDenseRestoreArray(EVAL,&Earray);
141: //VecView(eval,PETSC_VIEWER_STDOUT_WORLD);
142:
143: MatDiagonalScale(C2,NULL,eval); /* C2 = B*X*eval */
144: MatAXPY(C1,-1.0,C2,SAME_NONZERO_PATTERN); /* C1 = - C2 + C1 */
145: MatNorm(C1,NORM_FROBENIUS,&norm);
146: if (norm > 1.e-14) {
147: if (!rank) printf(" Warning: || A*X - B*X*We || = %g\n",norm);
148: }
150: MatDestroy(&C1);
151: MatDestroy(&C2);
152: VecDestroy(&eval);
153: MatDestroy(&EVAL);
154: MatDestroy(&X);
156: }
157: if (!isElemental) {
158: MatDestroy(&Ae);
159: MatDestroy(&Be);
161: /* Test MAT_REUSE_MATRIX which is only supported for inplace conversion */
162: MatConvert(A, MATELEMENTAL, MAT_REUSE_MATRIX, &A);
163: //MatView(A,PETSC_VIEWER_STDOUT_WORLD);
164: }
166: MatDestroy(&We);
167: MatDestroy(&Xe);
168: MatDestroy(&A);
169: MatDestroy(&B);
170: PetscFinalize();
171: return 0;
172: }