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