1 2 static char help[] = "Tests MatGetRowMax(), MatGetRowMin(), MatGetRowMaxAbs()\n"; 3 4 #include <petscmat.h> 5 6 #define M 5 7 #define N 6 8 9 int main(int argc,char **args) 10 { 11 Mat A; 12 Vec min,max,maxabs,e; 13 PetscInt m,n,j,imin[M],imax[M],imaxabs[M],indices[N],row,testcase=0; 14 PetscScalar values[N]; 15 PetscErrorCode ierr; 16 PetscMPIInt size,rank; 17 PetscReal enorm; 18 19 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 20 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 21 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 22 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-testcase",&testcase,NULL)); 23 24 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 25 if (testcase == 1) { /* proc[0] holds entire A and other processes have no entry */ 26 if (rank == 0) { 27 CHKERRQ(MatSetSizes(A,M,N,PETSC_DECIDE,PETSC_DECIDE)); 28 } else { 29 CHKERRQ(MatSetSizes(A,0,0,PETSC_DECIDE,PETSC_DECIDE)); 30 } 31 testcase = 0; 32 } else { 33 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N)); 34 } 35 CHKERRQ(MatSetFromOptions(A)); 36 CHKERRQ(MatSetUp(A)); 37 38 if (rank == 0) { /* proc[0] sets matrix A */ 39 for (j=0; j<N; j++) indices[j] = j; 40 switch (testcase) { 41 case 1: /* see testcast 0 */ 42 break; 43 case 2: 44 row = 0; 45 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -4.0; values[4] = 1.0; values[5] = 1.0; 46 CHKERRQ(MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES)); 47 row = 2; 48 indices[0] = 0; indices[1] = 3; indices[2] = 5; 49 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 50 CHKERRQ(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 51 row = 3; 52 indices[0] = 0; indices[1] = 1; indices[2] = 4; 53 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 54 CHKERRQ(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 55 row = 4; 56 indices[0] = 0; indices[1] = 1; indices[2] = 2; 57 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 58 CHKERRQ(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 59 break; 60 case 3: 61 row = 0; 62 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 63 CHKERRQ(MatSetValues(A,1,&row,3,indices+1,values,INSERT_VALUES)); 64 row = 1; 65 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 66 CHKERRQ(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 67 row = 2; 68 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 69 CHKERRQ(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 70 row = 3; 71 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -1.0; 72 CHKERRQ(MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES)); 73 row = 4; 74 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -1.0; 75 CHKERRQ(MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES)); 76 break; 77 78 default: 79 row = 0; 80 values[0] = -1.0; values[1] = 0.0; values[2] = 1.0; values[3] = 3.0; values[4] = 4.0; values[5] = -5.0; 81 CHKERRQ(MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES)); 82 row = 1; 83 CHKERRQ(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 84 row = 3; 85 CHKERRQ(MatSetValues(A,1,&row,1,indices+4,values+4,INSERT_VALUES)); 86 row = 4; 87 CHKERRQ(MatSetValues(A,1,&row,2,indices+4,values+4,INSERT_VALUES)); 88 } 89 } 90 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 91 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 92 CHKERRQ(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); 93 94 CHKERRQ(MatGetLocalSize(A, &m,&n)); 95 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&min)); 96 CHKERRQ(VecSetSizes(min,m,PETSC_DECIDE)); 97 CHKERRQ(VecSetFromOptions(min)); 98 CHKERRQ(VecDuplicate(min,&max)); 99 CHKERRQ(VecDuplicate(min,&maxabs)); 100 CHKERRQ(VecDuplicate(min,&e)); 101 102 /* MatGetRowMax() */ 103 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMax\n")); 104 CHKERRQ(MatGetRowMax(A,max,NULL)); 105 CHKERRQ(MatGetRowMax(A,max,imax)); 106 CHKERRQ(VecView(max,PETSC_VIEWER_STDOUT_WORLD)); 107 CHKERRQ(VecGetLocalSize(max,&n)); 108 CHKERRQ(PetscIntView(n,imax,PETSC_VIEWER_STDOUT_WORLD)); 109 110 /* MatGetRowMin() */ 111 CHKERRQ(MatScale(A,-1.0)); 112 CHKERRQ(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); 113 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMin\n")); 114 CHKERRQ(MatGetRowMin(A,min,NULL)); 115 CHKERRQ(MatGetRowMin(A,min,imin)); 116 117 CHKERRQ(VecWAXPY(e,1.0,max,min)); /* e = max + min */ 118 CHKERRQ(VecNorm(e,NORM_INFINITY,&enorm)); 119 PetscCheckFalse(enorm > PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"max+min > PETSC_MACHINE_EPSILON "); 120 for (j = 0; j < n; j++) { 121 PetscCheckFalse(imin[j] != imax[j],PETSC_COMM_SELF,PETSC_ERR_PLIB,"imin[%" PetscInt_FMT "] %" PetscInt_FMT " != imax %" PetscInt_FMT,j,imin[j],imax[j]); 122 } 123 124 /* MatGetRowMaxAbs() */ 125 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs\n")); 126 CHKERRQ(MatGetRowMaxAbs(A,maxabs,NULL)); 127 CHKERRQ(MatGetRowMaxAbs(A,maxabs,imaxabs)); 128 CHKERRQ(VecView(maxabs,PETSC_VIEWER_STDOUT_WORLD)); 129 CHKERRQ(PetscIntView(n,imaxabs,PETSC_VIEWER_STDOUT_WORLD)); 130 131 /* MatGetRowMinAbs() */ 132 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMinAbs\n")); 133 CHKERRQ(MatGetRowMinAbs(A,min,NULL)); 134 CHKERRQ(MatGetRowMinAbs(A,min,imin)); 135 CHKERRQ(VecView(min,PETSC_VIEWER_STDOUT_WORLD)); 136 CHKERRQ(PetscIntView(n,imin,PETSC_VIEWER_STDOUT_WORLD)); 137 138 if (size == 1) { 139 /* Test MatGetRowMax, MatGetRowMin and MatGetRowMaxAbs for SeqDense and MPIBAIJ matrix */ 140 Mat Adense; 141 Vec max_d,maxabs_d; 142 CHKERRQ(VecDuplicate(min,&max_d)); 143 CHKERRQ(VecDuplicate(min,&maxabs_d)); 144 145 CHKERRQ(MatScale(A,-1.0)); 146 CHKERRQ(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense)); 147 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMax for seqdense matrix\n")); 148 CHKERRQ(MatGetRowMax(Adense,max_d,imax)); 149 150 CHKERRQ(VecWAXPY(e,-1.0,max,max_d)); /* e = -max + max_d */ 151 CHKERRQ(VecNorm(e,NORM_INFINITY,&enorm)); 152 PetscCheckFalse(enorm > PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-max + max_d) %g > PETSC_MACHINE_EPSILON",(double)enorm); 153 154 CHKERRQ(MatScale(Adense,-1.0)); 155 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMin for seqdense matrix\n")); 156 CHKERRQ(MatGetRowMin(Adense,min,imin)); 157 158 CHKERRQ(VecWAXPY(e,1.0,max,min)); /* e = max + min */ 159 CHKERRQ(VecNorm(e,NORM_INFINITY,&enorm)); 160 PetscCheckFalse(enorm > PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"max+min > PETSC_MACHINE_EPSILON "); 161 for (j = 0; j < n; j++) { 162 if (imin[j] != imax[j]) { 163 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"imin[%" PetscInt_FMT "] %" PetscInt_FMT " != imax %" PetscInt_FMT " for seqdense matrix",j,imin[j],imax[j]); 164 } 165 } 166 167 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMaxAbs for seqdense matrix\n")); 168 CHKERRQ(MatGetRowMaxAbs(Adense,maxabs_d,imaxabs)); 169 CHKERRQ(VecWAXPY(e,-1.0,maxabs,maxabs_d)); /* e = -maxabs + maxabs_d */ 170 CHKERRQ(VecNorm(e,NORM_INFINITY,&enorm)); 171 PetscCheckFalse(enorm > PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON",(double)enorm); 172 173 CHKERRQ(MatDestroy(&Adense)); 174 CHKERRQ(VecDestroy(&max_d)); 175 CHKERRQ(VecDestroy(&maxabs_d)); 176 } 177 178 { /* BAIJ matrix */ 179 Mat B; 180 Vec maxabsB,maxabsB2; 181 PetscInt bs=2,*imaxabsB,*imaxabsB2,rstart,rend,cstart,cend,ncols,col,Brows[2],Bcols[2]; 182 const PetscInt *cols; 183 const PetscScalar *vals,*vals2; 184 PetscScalar Bvals[4]; 185 186 CHKERRQ(PetscMalloc2(M,&imaxabsB,bs*M,&imaxabsB2)); 187 188 /* bs = 1 */ 189 CHKERRQ(MatConvert(A,MATMPIBAIJ,MAT_INITIAL_MATRIX,&B)); 190 CHKERRQ(VecDuplicate(min,&maxabsB)); 191 CHKERRQ(MatGetRowMaxAbs(B,maxabsB,NULL)); 192 CHKERRQ(MatGetRowMaxAbs(B,maxabsB,imaxabsB)); 193 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs for BAIJ matrix\n")); 194 CHKERRQ(VecWAXPY(e,-1.0,maxabs,maxabsB)); /* e = -maxabs + maxabsB */ 195 CHKERRQ(VecNorm(e,NORM_INFINITY,&enorm)); 196 PetscCheckFalse(enorm > PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON",(double)enorm); 197 198 for (j = 0; j < n; j++) { 199 PetscCheckFalse(imaxabs[j] != imaxabsB[j],PETSC_COMM_SELF,PETSC_ERR_PLIB,"imaxabs[%" PetscInt_FMT "] %" PetscInt_FMT " != imaxabsB %" PetscInt_FMT,j,imin[j],imax[j]); 200 } 201 CHKERRQ(MatDestroy(&B)); 202 203 /* Test bs = 2: Create B with bs*bs block structure of A */ 204 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&maxabsB2)); 205 CHKERRQ(VecSetSizes(maxabsB2,bs*m,PETSC_DECIDE)); 206 CHKERRQ(VecSetFromOptions(maxabsB2)); 207 208 CHKERRQ(MatGetOwnershipRange(A,&rstart,&rend)); 209 CHKERRQ(MatGetOwnershipRangeColumn(A,&cstart,&cend)); 210 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&B)); 211 CHKERRQ(MatSetSizes(B,bs*(rend-rstart),bs*(cend-cstart),PETSC_DECIDE,PETSC_DECIDE)); 212 CHKERRQ(MatSetFromOptions(B)); 213 CHKERRQ(MatSetUp(B)); 214 215 for (row=rstart; row<rend; row++) { 216 CHKERRQ(MatGetRow(A,row,&ncols,&cols,&vals)); 217 for (col=0; col<ncols; col++) { 218 for (j=0; j<bs; j++) { 219 Brows[j] = bs*row + j; 220 Bcols[j] = bs*cols[col]+j; 221 } 222 for (j=0; j<bs*bs; j++) Bvals[j] = vals[col]; 223 CHKERRQ(MatSetValues(B,bs,Brows,bs,Bcols,Bvals,INSERT_VALUES)); 224 } 225 CHKERRQ(MatRestoreRow(A,row,&ncols,&cols,&vals)); 226 } 227 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 228 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 229 230 CHKERRQ(MatGetRowMaxAbs(B,maxabsB2,imaxabsB2)); 231 232 /* Check maxabsB2 and imaxabsB2 */ 233 CHKERRQ(VecGetArrayRead(maxabsB,&vals)); 234 CHKERRQ(VecGetArrayRead(maxabsB2,&vals2)); 235 for (row=0; row<m; row++) { 236 if (PetscAbsScalar(vals[row] - vals2[bs*row]) > PETSC_MACHINE_EPSILON) 237 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"row %" PetscInt_FMT " maxabsB != maxabsB2",row); 238 } 239 CHKERRQ(VecRestoreArrayRead(maxabsB,&vals)); 240 CHKERRQ(VecRestoreArrayRead(maxabsB2,&vals2)); 241 242 for (col=0; col<n; col++) { 243 if (imaxabsB[col] != imaxabsB2[bs*col]/bs) 244 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"col %" PetscInt_FMT " imaxabsB != imaxabsB2",col); 245 } 246 CHKERRQ(VecDestroy(&maxabsB)); 247 CHKERRQ(MatDestroy(&B)); 248 CHKERRQ(VecDestroy(&maxabsB2)); 249 CHKERRQ(PetscFree2(imaxabsB,imaxabsB2)); 250 } 251 252 CHKERRQ(VecDestroy(&min)); 253 CHKERRQ(VecDestroy(&max)); 254 CHKERRQ(VecDestroy(&maxabs)); 255 CHKERRQ(VecDestroy(&e)); 256 CHKERRQ(MatDestroy(&A)); 257 ierr = PetscFinalize(); 258 return ierr; 259 } 260 261 /*TEST 262 263 test: 264 output_file: output/ex114.out 265 266 test: 267 suffix: 2 268 args: -testcase 1 269 output_file: output/ex114.out 270 271 test: 272 suffix: 3 273 args: -testcase 2 274 output_file: output/ex114_3.out 275 276 test: 277 suffix: 4 278 args: -testcase 3 279 output_file: output/ex114_4.out 280 281 test: 282 suffix: 5 283 nsize: 3 284 args: -testcase 0 285 output_file: output/ex114_5.out 286 287 test: 288 suffix: 6 289 nsize: 3 290 args: -testcase 1 291 output_file: output/ex114_6.out 292 293 test: 294 suffix: 7 295 nsize: 3 296 args: -testcase 2 297 output_file: output/ex114_7.out 298 299 test: 300 suffix: 8 301 nsize: 3 302 args: -testcase 3 303 output_file: output/ex114_8.out 304 305 TEST*/ 306