Actual source code: ex174.cxx
petsc-3.7.3 2016-08-01
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: 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,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,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 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not supported yet");
102: }
103: MatConvert(A, MATELEMENTAL, MAT_INITIAL_MATRIX, &Ae);
104: MatConvert(B, MATELEMENTAL, MAT_INITIAL_MATRIX, &Be);
106: /* Test accuracy */
107: MatMultEqual(A,Ae,5,&flg);
108: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"A != A_elemental.");
109: MatMultEqual(B,Be,5,&flg);
110: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"B != B_elemental.");
111: }
113: /* Test MatElementalHermitianGenDefEig() */
114: if (!rank) printf(" Compute Ax = lambda Bx... \n");
115: vl = -0.8, vu = -0.7;
116: PetscOptionsGetReal(NULL,NULL,"-vl",&vl,NULL);
117: PetscOptionsGetReal(NULL,NULL,"-vu",&vu,NULL);
118: subset.rangeSubset = PETSC_TRUE;
119: subset.lowerBound = vl;
120: subset.upperBound = vu;
121: MatElementalHermitianGenDefEig(eigtype,uplo,Ae,Be,&We,&Xe,sort,subset,ctrl);
122: //MatView(We,PETSC_VIEWER_STDOUT_WORLD);
124: /* Check || A*X - B*X*We || */
125: if (isAij) {
126: if (!rank) printf(" Convert Elemental matrices We and Xe into MATDENSE matrices... \n");
127: MatConvert(We,MATDENSE,MAT_INITIAL_MATRIX,&EVAL); /* EVAL is a Mx1 matrix */
128: MatConvert(Xe,MATDENSE,MAT_INITIAL_MATRIX,&X);
129: //MatView(EVAL,PETSC_VIEWER_STDOUT_WORLD);
131: MatMatMult(A,X,MAT_INITIAL_MATRIX,1.0,&C1); /* C1 = A*X */
132: MatMatMult(B,X,MAT_INITIAL_MATRIX,1.0,&C2); /* C2 = B*X */
133:
134: /* Get vector eval from matrix EVAL for MatDiagonalScale() */
135: MatGetLocalSize(EVAL,&m,NULL);
136: MatDenseGetArray(EVAL,&Earray);
137: VecCreateMPIWithArray(PetscObjectComm((PetscObject)EVAL),1,m,PETSC_DECIDE,Earray,&eval);
138: MatDenseRestoreArray(EVAL,&Earray);
139: //VecView(eval,PETSC_VIEWER_STDOUT_WORLD);
140:
141: MatDiagonalScale(C2,NULL,eval); /* C2 = B*X*eval */
142: MatAXPY(C1,-1.0,C2,SAME_NONZERO_PATTERN); /* C1 = - C2 + C1 */
143: MatNorm(C1,NORM_FROBENIUS,&norm);
144: if (norm > 1.e-14) {
145: if (!rank) printf(" Warning: || A*X - B*X*We || = %g\n",norm);
146: }
148: MatDestroy(&C1);
149: MatDestroy(&C2);
150: VecDestroy(&eval);
151: MatDestroy(&EVAL);
152: MatDestroy(&X);
154: }
155: if (!isElemental) {
156: MatDestroy(&Ae);
157: MatDestroy(&Be);
159: /* Test MAT_REUSE_MATRIX which is only supported for inplace conversion */
160: MatConvert(A, MATELEMENTAL, MAT_INPLACE_MATRIX, &A);
161: //MatView(A,PETSC_VIEWER_STDOUT_WORLD);
162: }
164: MatDestroy(&We);
165: MatDestroy(&Xe);
166: MatDestroy(&A);
167: MatDestroy(&B);
168: PetscFinalize();
169: return 0;
170: }