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