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