1 2 static char help[] = "Tests the use of interface functions for MATIS matrices.\n\ 3 This example tests: MatZeroRows(), MatZeroRowsLocal(), MatView(), MatDuplicate(),\n\ 4 MatCopy(), MatCreateSubMatrix(), MatGetLocalSubMatrix(), MatAXPY(), MatShift()\n\ 5 MatDiagonalSet(), MatTranspose() and MatPtAP(). It also tests some\n\ 6 conversion routines.\n"; 7 8 #include <petscmat.h> 9 10 PetscErrorCode TestMatZeroRows(Mat,Mat,PetscBool,IS,PetscScalar); 11 PetscErrorCode CheckMat(Mat,Mat,PetscBool,const char*); 12 13 int main(int argc,char **args) 14 { 15 Mat A,B,A2,B2,T; 16 Mat Aee,Aeo,Aoe,Aoo; 17 Mat *mats; 18 Vec x,y; 19 MatInfo info; 20 ISLocalToGlobalMapping cmap,rmap; 21 IS is,is2,reven,rodd,ceven,codd; 22 IS *rows,*cols; 23 MatType lmtype; 24 PetscScalar diag = 2.; 25 PetscInt n,m,i,lm,ln; 26 PetscInt rst,ren,cst,cen,nr,nc; 27 PetscMPIInt rank,size; 28 PetscBool testT,squaretest,isaij; 29 PetscBool permute = PETSC_FALSE, negmap = PETSC_FALSE, repmap = PETSC_FALSE; 30 PetscBool diffmap = PETSC_TRUE, symmetric = PETSC_FALSE, issymmetric; 31 PetscErrorCode ierr; 32 33 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 34 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 35 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 36 m = n = 2*size; 37 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-symmetric",&symmetric,NULL)); 38 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 39 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 40 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-negmap",&negmap,NULL)); 41 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-repmap",&repmap,NULL)); 42 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-permmap",&permute,NULL)); 43 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-diffmap",&diffmap,NULL)); 44 PetscCheckFalse(size > 1 && m < 4,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of rows should be larger or equal 4 for parallel runs"); 45 PetscCheckFalse(size == 1 && m < 2,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of rows should be larger or equal 2 for uniprocessor runs"); 46 PetscCheckFalse(n < 2,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of cols should be larger or equal 2"); 47 48 /* create a MATIS matrix */ 49 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 50 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n)); 51 CHKERRQ(MatSetType(A,MATIS)); 52 CHKERRQ(MatSetFromOptions(A)); 53 if (!negmap && !repmap) { 54 /* This is not the proper setting for MATIS for finite elements, it is just used to test the routines 55 Here we use a one-to-one correspondence between local row/column spaces and global row/column spaces 56 Equivalent to passing NULL for the mapping */ 57 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,n,0,1,&is)); 58 } else if (negmap && !repmap) { /* non repeated but with negative indices */ 59 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,n+2,-2,1,&is)); 60 } else if (!negmap && repmap) { /* non negative but repeated indices */ 61 IS isl[2]; 62 63 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,n,0,1,&isl[0])); 64 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,n,n-1,-1,&isl[1])); 65 CHKERRQ(ISConcatenate(PETSC_COMM_WORLD,2,isl,&is)); 66 CHKERRQ(ISDestroy(&isl[0])); 67 CHKERRQ(ISDestroy(&isl[1])); 68 } else { /* negative and repeated indices */ 69 IS isl[2]; 70 71 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,n+1,-1,1,&isl[0])); 72 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,n+1,n-1,-1,&isl[1])); 73 CHKERRQ(ISConcatenate(PETSC_COMM_WORLD,2,isl,&is)); 74 CHKERRQ(ISDestroy(&isl[0])); 75 CHKERRQ(ISDestroy(&isl[1])); 76 } 77 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&cmap)); 78 CHKERRQ(ISDestroy(&is)); 79 80 if (m != n || diffmap) { 81 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,m,permute ? m-1 : 0,permute ? -1 : 1,&is)); 82 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&rmap)); 83 CHKERRQ(ISDestroy(&is)); 84 } else { 85 CHKERRQ(PetscObjectReference((PetscObject)cmap)); 86 rmap = cmap; 87 } 88 89 CHKERRQ(MatSetLocalToGlobalMapping(A,rmap,cmap)); 90 CHKERRQ(MatISStoreL2L(A,PETSC_FALSE)); 91 CHKERRQ(MatISSetPreallocation(A,3,NULL,3,NULL)); 92 CHKERRQ(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,(PetscBool)!(repmap || negmap))); /* I do not want to precompute the pattern */ 93 CHKERRQ(ISLocalToGlobalMappingGetSize(rmap,&lm)); 94 CHKERRQ(ISLocalToGlobalMappingGetSize(cmap,&ln)); 95 for (i=0; i<lm; i++) { 96 PetscScalar v[3]; 97 PetscInt cols[3]; 98 99 cols[0] = (i-1+n)%n; 100 cols[1] = i%n; 101 cols[2] = (i+1)%n; 102 v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1); 103 v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1); 104 v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1); 105 CHKERRQ(ISGlobalToLocalMappingApply(cmap,IS_GTOLM_MASK,3,cols,NULL,cols)); 106 CHKERRQ(MatSetValuesLocal(A,1,&i,3,cols,v,ADD_VALUES)); 107 } 108 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 109 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 110 111 /* activate tests for square matrices with same maps only */ 112 CHKERRQ(MatHasCongruentLayouts(A,&squaretest)); 113 if (squaretest && rmap != cmap) { 114 PetscInt nr, nc; 115 116 CHKERRQ(ISLocalToGlobalMappingGetSize(rmap,&nr)); 117 CHKERRQ(ISLocalToGlobalMappingGetSize(cmap,&nc)); 118 if (nr != nc) squaretest = PETSC_FALSE; 119 else { 120 const PetscInt *idxs1,*idxs2; 121 122 CHKERRQ(ISLocalToGlobalMappingGetIndices(rmap,&idxs1)); 123 CHKERRQ(ISLocalToGlobalMappingGetIndices(cmap,&idxs2)); 124 CHKERRQ(PetscArraycmp(idxs1,idxs2,nr,&squaretest)); 125 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(rmap,&idxs1)); 126 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(cmap,&idxs2)); 127 } 128 CHKERRMPI(MPIU_Allreduce(MPI_IN_PLACE,&squaretest,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A))); 129 } 130 131 /* test MatISGetLocalMat */ 132 CHKERRQ(MatISGetLocalMat(A,&B)); 133 CHKERRQ(MatGetType(B,&lmtype)); 134 135 /* test MatGetInfo */ 136 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatGetInfo\n")); 137 CHKERRQ(MatGetInfo(A,MAT_LOCAL,&info)); 138 CHKERRQ(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 139 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"Process %2d: %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",PetscGlobalRank,(PetscInt)info.nz_used, 140 (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);CHKERRQ(ierr); 141 CHKERRQ(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 142 CHKERRQ(MatGetInfo(A,MAT_GLOBAL_MAX,&info)); 143 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalMax : %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",(PetscInt)info.nz_used, 144 (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);CHKERRQ(ierr); 145 CHKERRQ(MatGetInfo(A,MAT_GLOBAL_SUM,&info)); 146 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalSum : %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",(PetscInt)info.nz_used, 147 (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);CHKERRQ(ierr); 148 149 /* test MatIsSymmetric */ 150 CHKERRQ(MatIsSymmetric(A,0.0,&issymmetric)); 151 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatIsSymmetric: %d\n",issymmetric)); 152 153 /* Create a MPIAIJ matrix, same as A */ 154 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&B)); 155 CHKERRQ(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n)); 156 CHKERRQ(MatSetType(B,MATAIJ)); 157 CHKERRQ(MatSetFromOptions(B)); 158 CHKERRQ(MatSetUp(B)); 159 CHKERRQ(MatSetLocalToGlobalMapping(B,rmap,cmap)); 160 CHKERRQ(MatMPIAIJSetPreallocation(B,3,NULL,3,NULL)); 161 CHKERRQ(MatMPIBAIJSetPreallocation(B,1,3,NULL,3,NULL)); 162 #if defined(PETSC_HAVE_HYPRE) 163 CHKERRQ(MatHYPRESetPreallocation(B,3,NULL,3,NULL)); 164 #endif 165 CHKERRQ(MatISSetPreallocation(B,3,NULL,3,NULL)); 166 CHKERRQ(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,(PetscBool)!(repmap || negmap))); /* I do not want to precompute the pattern */ 167 for (i=0; i<lm; i++) { 168 PetscScalar v[3]; 169 PetscInt cols[3]; 170 171 cols[0] = (i-1+n)%n; 172 cols[1] = i%n; 173 cols[2] = (i+1)%n; 174 v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1); 175 v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1); 176 v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1); 177 CHKERRQ(ISGlobalToLocalMappingApply(cmap,IS_GTOLM_MASK,3,cols,NULL,cols)); 178 CHKERRQ(MatSetValuesLocal(B,1,&i,3,cols,v,ADD_VALUES)); 179 } 180 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 181 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 182 183 /* test MatView */ 184 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatView\n")); 185 CHKERRQ(MatView(A,NULL)); 186 CHKERRQ(MatView(B,NULL)); 187 188 /* test CheckMat */ 189 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test CheckMat\n")); 190 CHKERRQ(CheckMat(A,B,PETSC_FALSE,"CheckMat")); 191 192 /* test MatDuplicate and MatAXPY */ 193 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatDuplicate and MatAXPY\n")); 194 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 195 CHKERRQ(CheckMat(A,A2,PETSC_FALSE,"MatDuplicate and MatAXPY")); 196 197 /* test MatConvert */ 198 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_IS_XAIJ\n")); 199 CHKERRQ(MatConvert(A2,MATAIJ,MAT_INITIAL_MATRIX,&B2)); 200 CHKERRQ(CheckMat(B,B2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_INITIAL_MATRIX")); 201 CHKERRQ(MatConvert(A2,MATAIJ,MAT_REUSE_MATRIX,&B2)); 202 CHKERRQ(CheckMat(B,B2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_REUSE_MATRIX")); 203 CHKERRQ(MatConvert(A2,MATAIJ,MAT_INPLACE_MATRIX,&A2)); 204 CHKERRQ(CheckMat(B,A2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_INPLACE_MATRIX")); 205 CHKERRQ(MatDestroy(&A2)); 206 CHKERRQ(MatDestroy(&B2)); 207 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_XAIJ_IS\n")); 208 CHKERRQ(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 209 CHKERRQ(MatConvert(B2,MATIS,MAT_INITIAL_MATRIX,&A2)); 210 CHKERRQ(CheckMat(A,A2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_INITIAL_MATRIX")); 211 CHKERRQ(MatConvert(B2,MATIS,MAT_REUSE_MATRIX,&A2)); 212 CHKERRQ(CheckMat(A,A2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_REUSE_MATRIX")); 213 CHKERRQ(MatConvert(B2,MATIS,MAT_INPLACE_MATRIX,&B2)); 214 CHKERRQ(CheckMat(A,B2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_INPLACE_MATRIX")); 215 CHKERRQ(MatDestroy(&A2)); 216 CHKERRQ(MatDestroy(&B2)); 217 CHKERRQ(PetscStrcmp(lmtype,MATSEQAIJ,&isaij)); 218 if (size == 1 && isaij) { /* tests special code paths in MatConvert_IS_XAIJ */ 219 PetscInt ri, ci, rr[3] = {0,1,0}, cr[4] = {1,2,0,1}, rk[3] = {0,2,1}, ck[4] = {1,0,3,2}; 220 ISLocalToGlobalMapping tcmap,trmap; 221 222 for (ri = 0; ri < 2; ri++) { 223 PetscInt *r; 224 225 r = (PetscInt*)(ri == 0 ? rr : rk); 226 for (ci = 0; ci < 2; ci++) { 227 PetscInt *c,rb,cb; 228 229 c = (PetscInt*)(ci == 0 ? cr : ck); 230 for (rb = 1; rb < 4; rb++) { 231 CHKERRQ(ISCreateBlock(PETSC_COMM_SELF,rb,3,r,PETSC_COPY_VALUES,&is)); 232 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&trmap)); 233 CHKERRQ(ISDestroy(&is)); 234 for (cb = 1; cb < 4; cb++) { 235 Mat T,lT,T2; 236 char testname[256]; 237 238 CHKERRQ(PetscSNPrintf(testname,sizeof(testname),"MatConvert_IS_XAIJ special case (%" PetscInt_FMT " %" PetscInt_FMT ", bs %" PetscInt_FMT " %" PetscInt_FMT ")",ri,ci,rb,cb)); 239 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test %s\n",testname)); 240 241 CHKERRQ(ISCreateBlock(PETSC_COMM_SELF,cb,4,c,PETSC_COPY_VALUES,&is)); 242 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&tcmap)); 243 CHKERRQ(ISDestroy(&is)); 244 245 CHKERRQ(MatCreate(PETSC_COMM_SELF,&T)); 246 CHKERRQ(MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,rb*3,cb*4)); 247 CHKERRQ(MatSetType(T,MATIS)); 248 CHKERRQ(MatSetLocalToGlobalMapping(T,trmap,tcmap)); 249 CHKERRQ(ISLocalToGlobalMappingDestroy(&tcmap)); 250 CHKERRQ(MatISGetLocalMat(T,&lT)); 251 CHKERRQ(MatSetType(lT,MATSEQAIJ)); 252 CHKERRQ(MatSeqAIJSetPreallocation(lT,cb*4,NULL)); 253 CHKERRQ(MatSetRandom(lT,NULL)); 254 CHKERRQ(MatConvert(lT,lmtype,MAT_INPLACE_MATRIX,&lT)); 255 CHKERRQ(MatISRestoreLocalMat(T,&lT)); 256 CHKERRQ(MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY)); 257 CHKERRQ(MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY)); 258 259 CHKERRQ(MatConvert(T,MATAIJ,MAT_INITIAL_MATRIX,&T2)); 260 CHKERRQ(CheckMat(T,T2,PETSC_TRUE,"MAT_INITIAL_MATRIX")); 261 CHKERRQ(MatConvert(T,MATAIJ,MAT_REUSE_MATRIX,&T2)); 262 CHKERRQ(CheckMat(T,T2,PETSC_TRUE,"MAT_REUSE_MATRIX")); 263 CHKERRQ(MatDestroy(&T2)); 264 CHKERRQ(MatDuplicate(T,MAT_COPY_VALUES,&T2)); 265 CHKERRQ(MatConvert(T2,MATAIJ,MAT_INPLACE_MATRIX,&T2)); 266 CHKERRQ(CheckMat(T,T2,PETSC_TRUE,"MAT_INPLACE_MATRIX")); 267 CHKERRQ(MatDestroy(&T)); 268 CHKERRQ(MatDestroy(&T2)); 269 } 270 CHKERRQ(ISLocalToGlobalMappingDestroy(&trmap)); 271 } 272 } 273 } 274 } 275 276 /* test MatDiagonalScale */ 277 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalScale\n")); 278 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 279 CHKERRQ(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 280 CHKERRQ(MatCreateVecs(A,&x,&y)); 281 CHKERRQ(VecSetRandom(x,NULL)); 282 if (issymmetric) { 283 CHKERRQ(VecCopy(x,y)); 284 } else { 285 CHKERRQ(VecSetRandom(y,NULL)); 286 CHKERRQ(VecScale(y,8.)); 287 } 288 CHKERRQ(MatDiagonalScale(A2,y,x)); 289 CHKERRQ(MatDiagonalScale(B2,y,x)); 290 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalScale")); 291 CHKERRQ(MatDestroy(&A2)); 292 CHKERRQ(MatDestroy(&B2)); 293 CHKERRQ(VecDestroy(&x)); 294 CHKERRQ(VecDestroy(&y)); 295 296 /* test MatPtAP (A IS and B AIJ) */ 297 if (isaij && m == n) { 298 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatPtAP\n")); 299 CHKERRQ(MatISStoreL2L(A,PETSC_TRUE)); 300 CHKERRQ(MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&A2)); 301 CHKERRQ(MatPtAP(B,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B2)); 302 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"MatPtAP MAT_INITIAL_MATRIX")); 303 CHKERRQ(MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&A2)); 304 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"MatPtAP MAT_REUSE_MATRIX")); 305 CHKERRQ(MatDestroy(&A2)); 306 CHKERRQ(MatDestroy(&B2)); 307 } 308 309 /* test MatGetLocalSubMatrix */ 310 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatGetLocalSubMatrix\n")); 311 CHKERRQ(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&A2)); 312 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,lm/2+lm%2,0,2,&reven)); 313 CHKERRQ(ISComplement(reven,0,lm,&rodd)); 314 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,ln/2+ln%2,0,2,&ceven)); 315 CHKERRQ(ISComplement(ceven,0,ln,&codd)); 316 CHKERRQ(MatGetLocalSubMatrix(A2,reven,ceven,&Aee)); 317 CHKERRQ(MatGetLocalSubMatrix(A2,reven,codd,&Aeo)); 318 CHKERRQ(MatGetLocalSubMatrix(A2,rodd,ceven,&Aoe)); 319 CHKERRQ(MatGetLocalSubMatrix(A2,rodd,codd,&Aoo)); 320 for (i=0; i<lm; i++) { 321 PetscInt j,je,jo,colse[3], colso[3]; 322 PetscScalar ve[3], vo[3]; 323 PetscScalar v[3]; 324 PetscInt cols[3]; 325 PetscInt row = i/2; 326 327 cols[0] = (i-1+n)%n; 328 cols[1] = i%n; 329 cols[2] = (i+1)%n; 330 v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1); 331 v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1); 332 v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1); 333 CHKERRQ(ISGlobalToLocalMappingApply(cmap,IS_GTOLM_MASK,3,cols,NULL,cols)); 334 for (j=0,je=0,jo=0;j<3;j++) { 335 if (cols[j]%2) { 336 vo[jo] = v[j]; 337 colso[jo++] = cols[j]/2; 338 } else { 339 ve[je] = v[j]; 340 colse[je++] = cols[j]/2; 341 } 342 } 343 if (i%2) { 344 CHKERRQ(MatSetValuesLocal(Aoe,1,&row,je,colse,ve,ADD_VALUES)); 345 CHKERRQ(MatSetValuesLocal(Aoo,1,&row,jo,colso,vo,ADD_VALUES)); 346 } else { 347 CHKERRQ(MatSetValuesLocal(Aee,1,&row,je,colse,ve,ADD_VALUES)); 348 CHKERRQ(MatSetValuesLocal(Aeo,1,&row,jo,colso,vo,ADD_VALUES)); 349 } 350 } 351 CHKERRQ(MatRestoreLocalSubMatrix(A2,reven,ceven,&Aee)); 352 CHKERRQ(MatRestoreLocalSubMatrix(A2,reven,codd,&Aeo)); 353 CHKERRQ(MatRestoreLocalSubMatrix(A2,rodd,ceven,&Aoe)); 354 CHKERRQ(MatRestoreLocalSubMatrix(A2,rodd,codd,&Aoo)); 355 CHKERRQ(ISDestroy(&reven)); 356 CHKERRQ(ISDestroy(&ceven)); 357 CHKERRQ(ISDestroy(&rodd)); 358 CHKERRQ(ISDestroy(&codd)); 359 CHKERRQ(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 360 CHKERRQ(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 361 CHKERRQ(MatAXPY(A2,-1.,A,SAME_NONZERO_PATTERN)); 362 CHKERRQ(CheckMat(A2,NULL,PETSC_FALSE,"MatGetLocalSubMatrix")); 363 CHKERRQ(MatDestroy(&A2)); 364 365 /* test MatConvert_Nest_IS */ 366 testT = PETSC_FALSE; 367 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_trans",&testT,NULL)); 368 369 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_Nest_IS\n")); 370 nr = 2; 371 nc = 2; 372 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-nr",&nr,NULL)); 373 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL)); 374 if (testT) { 375 CHKERRQ(MatGetOwnershipRange(A,&cst,&cen)); 376 CHKERRQ(MatGetOwnershipRangeColumn(A,&rst,&ren)); 377 } else { 378 CHKERRQ(MatGetOwnershipRange(A,&rst,&ren)); 379 CHKERRQ(MatGetOwnershipRangeColumn(A,&cst,&cen)); 380 } 381 CHKERRQ(PetscMalloc3(nr,&rows,nc,&cols,2*nr*nc,&mats)); 382 for (i=0;i<nr*nc;i++) { 383 if (testT) { 384 CHKERRQ(MatCreateTranspose(A,&mats[i])); 385 CHKERRQ(MatTranspose(B,MAT_INITIAL_MATRIX,&mats[i+nr*nc])); 386 } else { 387 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&mats[i])); 388 CHKERRQ(MatDuplicate(B,MAT_COPY_VALUES,&mats[i+nr*nc])); 389 } 390 } 391 for (i=0;i<nr;i++) { 392 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,ren-rst,i+rst,nr,&rows[i])); 393 } 394 for (i=0;i<nc;i++) { 395 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,cen-cst,i+cst,nc,&cols[i])); 396 } 397 CHKERRQ(MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats,&A2)); 398 CHKERRQ(MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats+nr*nc,&B2)); 399 for (i=0;i<nr;i++) { 400 CHKERRQ(ISDestroy(&rows[i])); 401 } 402 for (i=0;i<nc;i++) { 403 CHKERRQ(ISDestroy(&cols[i])); 404 } 405 for (i=0;i<2*nr*nc;i++) { 406 CHKERRQ(MatDestroy(&mats[i])); 407 } 408 CHKERRQ(PetscFree3(rows,cols,mats)); 409 CHKERRQ(MatConvert(B2,MATAIJ,MAT_INITIAL_MATRIX,&T)); 410 CHKERRQ(MatDestroy(&B2)); 411 CHKERRQ(MatConvert(A2,MATIS,MAT_INITIAL_MATRIX,&B2)); 412 CHKERRQ(CheckMat(B2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_INITIAL_MATRIX")); 413 CHKERRQ(MatConvert(A2,MATIS,MAT_REUSE_MATRIX,&B2)); 414 CHKERRQ(CheckMat(B2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_REUSE_MATRIX")); 415 CHKERRQ(MatDestroy(&B2)); 416 CHKERRQ(MatConvert(A2,MATIS,MAT_INPLACE_MATRIX,&A2)); 417 CHKERRQ(CheckMat(A2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_INPLACE_MATRIX")); 418 CHKERRQ(MatDestroy(&T)); 419 CHKERRQ(MatDestroy(&A2)); 420 421 /* test MatCreateSubMatrix */ 422 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatCreateSubMatrix\n")); 423 if (rank == 0) { 424 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,1,1,1,&is)); 425 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,2,0,1,&is2)); 426 } else if (rank == 1) { 427 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is)); 428 if (n > 3) { 429 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,1,3,1,&is2)); 430 } else { 431 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2)); 432 } 433 } else if (rank == 2 && n > 4) { 434 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is)); 435 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,n-4,4,1,&is2)); 436 } else { 437 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is)); 438 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2)); 439 } 440 CHKERRQ(MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&A2)); 441 CHKERRQ(MatCreateSubMatrix(B,is,is,MAT_INITIAL_MATRIX,&B2)); 442 CHKERRQ(CheckMat(A2,B2,PETSC_TRUE,"first MatCreateSubMatrix")); 443 444 CHKERRQ(MatCreateSubMatrix(A,is,is,MAT_REUSE_MATRIX,&A2)); 445 CHKERRQ(MatCreateSubMatrix(B,is,is,MAT_REUSE_MATRIX,&B2)); 446 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"reuse MatCreateSubMatrix")); 447 CHKERRQ(MatDestroy(&A2)); 448 CHKERRQ(MatDestroy(&B2)); 449 450 if (!issymmetric) { 451 CHKERRQ(MatCreateSubMatrix(A,is,is2,MAT_INITIAL_MATRIX,&A2)); 452 CHKERRQ(MatCreateSubMatrix(B,is,is2,MAT_INITIAL_MATRIX,&B2)); 453 CHKERRQ(MatCreateSubMatrix(A,is,is2,MAT_REUSE_MATRIX,&A2)); 454 CHKERRQ(MatCreateSubMatrix(B,is,is2,MAT_REUSE_MATRIX,&B2)); 455 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"second MatCreateSubMatrix")); 456 } 457 458 CHKERRQ(MatDestroy(&A2)); 459 CHKERRQ(MatDestroy(&B2)); 460 CHKERRQ(ISDestroy(&is)); 461 CHKERRQ(ISDestroy(&is2)); 462 463 /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */ 464 if (size > 1) { 465 if (rank == 0) { 466 PetscInt st,len; 467 468 st = (m+1)/2; 469 len = PetscMin(m/2,PetscMax(m-(m+1)/2-1,0)); 470 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,len,st,1,&is)); 471 } else { 472 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is)); 473 } 474 } else { 475 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is)); 476 } 477 478 if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */ 479 /* test MatDiagonalSet */ 480 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalSet\n")); 481 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 482 CHKERRQ(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 483 CHKERRQ(MatCreateVecs(A,NULL,&x)); 484 CHKERRQ(VecSetRandom(x,NULL)); 485 CHKERRQ(MatDiagonalSet(A2,x,INSERT_VALUES)); 486 CHKERRQ(MatDiagonalSet(B2,x,INSERT_VALUES)); 487 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalSet")); 488 CHKERRQ(VecDestroy(&x)); 489 CHKERRQ(MatDestroy(&A2)); 490 CHKERRQ(MatDestroy(&B2)); 491 492 /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */ 493 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatShift\n")); 494 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 495 CHKERRQ(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 496 CHKERRQ(MatShift(A2,2.0)); 497 CHKERRQ(MatShift(B2,2.0)); 498 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"MatShift")); 499 CHKERRQ(MatDestroy(&A2)); 500 CHKERRQ(MatDestroy(&B2)); 501 502 /* nonzero diag value is supported for square matrices only */ 503 CHKERRQ(TestMatZeroRows(A,B,PETSC_TRUE,is,diag)); 504 } 505 CHKERRQ(TestMatZeroRows(A,B,squaretest,is,0.0)); 506 CHKERRQ(ISDestroy(&is)); 507 508 /* test MatTranspose */ 509 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatTranspose\n")); 510 CHKERRQ(MatTranspose(A,MAT_INITIAL_MATRIX,&A2)); 511 CHKERRQ(MatTranspose(B,MAT_INITIAL_MATRIX,&B2)); 512 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"initial matrix MatTranspose")); 513 514 CHKERRQ(MatTranspose(A,MAT_REUSE_MATRIX,&A2)); 515 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (not in place) MatTranspose")); 516 CHKERRQ(MatDestroy(&A2)); 517 518 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 519 CHKERRQ(MatTranspose(A2,MAT_INPLACE_MATRIX,&A2)); 520 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (in place) MatTranspose")); 521 CHKERRQ(MatDestroy(&A2)); 522 523 CHKERRQ(MatTranspose(A,MAT_INITIAL_MATRIX,&A2)); 524 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (different type) MatTranspose")); 525 CHKERRQ(MatDestroy(&A2)); 526 CHKERRQ(MatDestroy(&B2)); 527 528 /* test MatISFixLocalEmpty */ 529 if (isaij) { 530 PetscInt r[2]; 531 532 r[0] = 0; 533 r[1] = PetscMin(m,n)-1; 534 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatISFixLocalEmpty\n")); 535 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 536 537 CHKERRQ(MatISFixLocalEmpty(A2,PETSC_TRUE)); 538 CHKERRQ(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 539 CHKERRQ(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 540 CHKERRQ(CheckMat(A2,B,PETSC_FALSE,"MatISFixLocalEmpty (null)")); 541 542 CHKERRQ(MatZeroRows(A2,2,r,0.0,NULL,NULL)); 543 CHKERRQ(MatViewFromOptions(A2,NULL,"-fixempty_view")); 544 CHKERRQ(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 545 CHKERRQ(MatZeroRows(B2,2,r,0.0,NULL,NULL)); 546 CHKERRQ(MatISFixLocalEmpty(A2,PETSC_TRUE)); 547 CHKERRQ(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 548 CHKERRQ(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 549 CHKERRQ(MatViewFromOptions(A2,NULL,"-fixempty_view")); 550 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows)")); 551 CHKERRQ(MatDestroy(&A2)); 552 553 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 554 CHKERRQ(MatZeroRows(A2,2,r,0.0,NULL,NULL)); 555 CHKERRQ(MatTranspose(A2,MAT_INPLACE_MATRIX,&A2)); 556 CHKERRQ(MatTranspose(B2,MAT_INPLACE_MATRIX,&B2)); 557 CHKERRQ(MatViewFromOptions(A2,NULL,"-fixempty_view")); 558 CHKERRQ(MatISFixLocalEmpty(A2,PETSC_TRUE)); 559 CHKERRQ(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 560 CHKERRQ(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 561 CHKERRQ(MatViewFromOptions(A2,NULL,"-fixempty_view")); 562 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (cols)")); 563 564 CHKERRQ(MatDestroy(&A2)); 565 CHKERRQ(MatDestroy(&B2)); 566 567 if (squaretest) { 568 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 569 CHKERRQ(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 570 CHKERRQ(MatZeroRowsColumns(A2,2,r,0.0,NULL,NULL)); 571 CHKERRQ(MatZeroRowsColumns(B2,2,r,0.0,NULL,NULL)); 572 CHKERRQ(MatViewFromOptions(A2,NULL,"-fixempty_view")); 573 CHKERRQ(MatISFixLocalEmpty(A2,PETSC_TRUE)); 574 CHKERRQ(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 575 CHKERRQ(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 576 CHKERRQ(MatViewFromOptions(A2,NULL,"-fixempty_view")); 577 CHKERRQ(CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows+cols)")); 578 CHKERRQ(MatDestroy(&A2)); 579 CHKERRQ(MatDestroy(&B2)); 580 } 581 582 } 583 584 /* test MatInvertBlockDiagonal 585 special cases for block-diagonal matrices */ 586 if (m == n) { 587 ISLocalToGlobalMapping map; 588 Mat Abd,Bbd; 589 IS is,bis; 590 const PetscScalar *isbd,*aijbd; 591 PetscScalar *vals; 592 const PetscInt *sts,*idxs; 593 PetscInt *idxs2,diff,perm,nl,bs,st,en,in; 594 PetscBool ok; 595 596 for (diff = 0; diff < 3; diff++) { 597 for (perm = 0; perm < 3; perm++) { 598 for (bs = 1; bs < 4; bs++) { 599 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatInvertBlockDiagonal blockdiag %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",n,diff,perm,bs)); 600 CHKERRQ(PetscMalloc1(bs*bs,&vals)); 601 CHKERRQ(MatGetOwnershipRanges(A,&sts)); 602 switch (diff) { 603 case 1: /* inverted layout by processes */ 604 in = 1; 605 st = sts[size - rank - 1]; 606 en = sts[size - rank]; 607 nl = en - st; 608 break; 609 case 2: /* round-robin layout */ 610 in = size; 611 st = rank; 612 nl = n/size; 613 if (rank < n%size) nl++; 614 break; 615 default: /* same layout */ 616 in = 1; 617 st = sts[rank]; 618 en = sts[rank + 1]; 619 nl = en - st; 620 break; 621 } 622 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,nl,st,in,&is)); 623 CHKERRQ(ISGetLocalSize(is,&nl)); 624 CHKERRQ(ISGetIndices(is,&idxs)); 625 CHKERRQ(PetscMalloc1(nl,&idxs2)); 626 for (i=0;i<nl;i++) { 627 switch (perm) { /* invert some of the indices */ 628 case 2: 629 idxs2[i] = rank%2 ? idxs[i] : idxs[nl-i-1]; 630 break; 631 case 1: 632 idxs2[i] = rank%2 ? idxs[nl-i-1] : idxs[i]; 633 break; 634 default: 635 idxs2[i] = idxs[i]; 636 break; 637 } 638 } 639 CHKERRQ(ISRestoreIndices(is,&idxs)); 640 CHKERRQ(ISCreateBlock(PETSC_COMM_WORLD,bs,nl,idxs2,PETSC_OWN_POINTER,&bis)); 641 CHKERRQ(ISLocalToGlobalMappingCreateIS(bis,&map)); 642 CHKERRQ(MatCreateIS(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,bs*n,bs*n,map,map,&Abd)); 643 CHKERRQ(ISLocalToGlobalMappingDestroy(&map)); 644 CHKERRQ(MatISSetPreallocation(Abd,bs,NULL,0,NULL)); 645 for (i=0;i<nl;i++) { 646 PetscInt b1,b2; 647 648 for (b1=0;b1<bs;b1++) { 649 for (b2=0;b2<bs;b2++) { 650 vals[b1*bs + b2] = i*bs*bs + b1*bs + b2 + 1 + (b1 == b2 ? 1.0 : 0); 651 } 652 } 653 CHKERRQ(MatSetValuesBlockedLocal(Abd,1,&i,1,&i,vals,INSERT_VALUES)); 654 } 655 CHKERRQ(MatAssemblyBegin(Abd,MAT_FINAL_ASSEMBLY)); 656 CHKERRQ(MatAssemblyEnd(Abd,MAT_FINAL_ASSEMBLY)); 657 CHKERRQ(MatConvert(Abd,MATAIJ,MAT_INITIAL_MATRIX,&Bbd)); 658 CHKERRQ(MatInvertBlockDiagonal(Abd,&isbd)); 659 CHKERRQ(MatInvertBlockDiagonal(Bbd,&aijbd)); 660 CHKERRQ(MatGetLocalSize(Bbd,&nl,NULL)); 661 ok = PETSC_TRUE; 662 for (i=0;i<nl/bs;i++) { 663 PetscInt b1,b2; 664 665 for (b1=0;b1<bs;b1++) { 666 for (b2=0;b2<bs;b2++) { 667 if (PetscAbsScalar(isbd[i*bs*bs+b1*bs + b2]-aijbd[i*bs*bs+b1*bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE; 668 if (!ok) { 669 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"[%d] ERROR block %" PetscInt_FMT ", entry %" PetscInt_FMT " %" PetscInt_FMT ": %g %g\n",rank,i,b1,b2,(double)PetscAbsScalar(isbd[i*bs*bs+b1*bs + b2]),(double)PetscAbsScalar(aijbd[i*bs*bs+b1*bs + b2]))); 670 break; 671 } 672 } 673 if (!ok) break; 674 } 675 if (!ok) break; 676 } 677 CHKERRQ(MatDestroy(&Abd)); 678 CHKERRQ(MatDestroy(&Bbd)); 679 CHKERRQ(PetscFree(vals)); 680 CHKERRQ(ISDestroy(&is)); 681 CHKERRQ(ISDestroy(&bis)); 682 } 683 } 684 } 685 } 686 /* free testing matrices */ 687 CHKERRQ(ISLocalToGlobalMappingDestroy(&cmap)); 688 CHKERRQ(ISLocalToGlobalMappingDestroy(&rmap)); 689 CHKERRQ(MatDestroy(&A)); 690 CHKERRQ(MatDestroy(&B)); 691 ierr = PetscFinalize(); 692 return ierr; 693 } 694 695 PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char* func) 696 { 697 Mat Bcheck; 698 PetscReal error; 699 700 PetscFunctionBeginUser; 701 if (!usemult && B) { 702 PetscBool hasnorm; 703 704 CHKERRQ(MatHasOperation(B,MATOP_NORM,&hasnorm)); 705 if (!hasnorm) usemult = PETSC_TRUE; 706 } 707 if (!usemult) { 708 if (B) { 709 MatType Btype; 710 711 CHKERRQ(MatGetType(B,&Btype)); 712 CHKERRQ(MatConvert(A,Btype,MAT_INITIAL_MATRIX,&Bcheck)); 713 } else { 714 CHKERRQ(MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck)); 715 } 716 if (B) { /* if B is present, subtract it */ 717 CHKERRQ(MatAXPY(Bcheck,-1.,B,DIFFERENT_NONZERO_PATTERN)); 718 } 719 CHKERRQ(MatNorm(Bcheck,NORM_INFINITY,&error)); 720 if (error > PETSC_SQRT_MACHINE_EPSILON) { 721 ISLocalToGlobalMapping rl2g,cl2g; 722 723 CHKERRQ(PetscObjectSetName((PetscObject)Bcheck,"Assembled Bcheck")); 724 CHKERRQ(MatView(Bcheck,NULL)); 725 if (B) { 726 CHKERRQ(PetscObjectSetName((PetscObject)B,"Assembled AIJ")); 727 CHKERRQ(MatView(B,NULL)); 728 CHKERRQ(MatDestroy(&Bcheck)); 729 CHKERRQ(MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck)); 730 CHKERRQ(PetscObjectSetName((PetscObject)Bcheck,"Assembled IS")); 731 CHKERRQ(MatView(Bcheck,NULL)); 732 } 733 CHKERRQ(MatDestroy(&Bcheck)); 734 CHKERRQ(PetscObjectSetName((PetscObject)A,"MatIS")); 735 CHKERRQ(MatView(A,NULL)); 736 CHKERRQ(MatGetLocalToGlobalMapping(A,&rl2g,&cl2g)); 737 CHKERRQ(ISLocalToGlobalMappingView(rl2g,NULL)); 738 CHKERRQ(ISLocalToGlobalMappingView(cl2g,NULL)); 739 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: %g",func,(double)error); 740 } 741 CHKERRQ(MatDestroy(&Bcheck)); 742 } else { 743 PetscBool ok,okt; 744 745 CHKERRQ(MatMultEqual(A,B,3,&ok)); 746 CHKERRQ(MatMultTransposeEqual(A,B,3,&okt)); 747 PetscCheckFalse(!ok || !okt,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: mult ok ? %d, multtranspose ok ? %d",func,ok,okt); 748 } 749 PetscFunctionReturn(0); 750 } 751 752 PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag) 753 { 754 Mat B,Bcheck,B2 = NULL,lB; 755 Vec x = NULL, b = NULL, b2 = NULL; 756 ISLocalToGlobalMapping l2gr,l2gc; 757 PetscReal error; 758 char diagstr[16]; 759 const PetscInt *idxs; 760 PetscInt rst,ren,i,n,N,d; 761 PetscMPIInt rank; 762 PetscBool miss,haszerorows; 763 764 PetscFunctionBeginUser; 765 if (diag == 0.) { 766 CHKERRQ(PetscStrcpy(diagstr,"zero")); 767 } else { 768 CHKERRQ(PetscStrcpy(diagstr,"nonzero")); 769 } 770 CHKERRQ(ISView(is,NULL)); 771 CHKERRQ(MatGetLocalToGlobalMapping(A,&l2gr,&l2gc)); 772 /* tests MatDuplicate and MatCopy */ 773 if (diag == 0.) { 774 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&B)); 775 } else { 776 CHKERRQ(MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B)); 777 CHKERRQ(MatCopy(A,B,SAME_NONZERO_PATTERN)); 778 } 779 CHKERRQ(MatISGetLocalMat(B,&lB)); 780 CHKERRQ(MatHasOperation(lB,MATOP_ZERO_ROWS,&haszerorows)); 781 if (squaretest && haszerorows) { 782 783 CHKERRQ(MatCreateVecs(B,&x,&b)); 784 CHKERRQ(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 785 CHKERRQ(VecSetLocalToGlobalMapping(b,l2gr)); 786 CHKERRQ(VecSetLocalToGlobalMapping(x,l2gc)); 787 CHKERRQ(VecSetRandom(x,NULL)); 788 CHKERRQ(VecSetRandom(b,NULL)); 789 /* mimic b[is] = x[is] */ 790 CHKERRQ(VecDuplicate(b,&b2)); 791 CHKERRQ(VecSetLocalToGlobalMapping(b2,l2gr)); 792 CHKERRQ(VecCopy(b,b2)); 793 CHKERRQ(ISGetLocalSize(is,&n)); 794 CHKERRQ(ISGetIndices(is,&idxs)); 795 CHKERRQ(VecGetSize(x,&N)); 796 for (i=0;i<n;i++) { 797 if (0 <= idxs[i] && idxs[i] < N) { 798 CHKERRQ(VecSetValue(b2,idxs[i],diag,INSERT_VALUES)); 799 CHKERRQ(VecSetValue(x,idxs[i],1.,INSERT_VALUES)); 800 } 801 } 802 CHKERRQ(VecAssemblyBegin(b2)); 803 CHKERRQ(VecAssemblyEnd(b2)); 804 CHKERRQ(VecAssemblyBegin(x)); 805 CHKERRQ(VecAssemblyEnd(x)); 806 CHKERRQ(ISRestoreIndices(is,&idxs)); 807 /* test ZeroRows on MATIS */ 808 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr)); 809 CHKERRQ(MatZeroRowsIS(B,is,diag,x,b)); 810 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRowsColumns (diag %s)\n",diagstr)); 811 CHKERRQ(MatZeroRowsColumnsIS(B2,is,diag,NULL,NULL)); 812 } else if (haszerorows) { 813 /* test ZeroRows on MATIS */ 814 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr)); 815 CHKERRQ(MatZeroRowsIS(B,is,diag,NULL,NULL)); 816 b = b2 = x = NULL; 817 } else { 818 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Skipping MatZeroRows (diag %s)\n",diagstr)); 819 b = b2 = x = NULL; 820 } 821 822 if (squaretest && haszerorows) { 823 CHKERRQ(VecAXPY(b2,-1.,b)); 824 CHKERRQ(VecNorm(b2,NORM_INFINITY,&error)); 825 PetscCheckFalse(error > PETSC_SQRT_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR IN ZEROROWS ON B %g (diag %s)",(double)error,diagstr); 826 } 827 828 /* test MatMissingDiagonal */ 829 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Test MatMissingDiagonal\n")); 830 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 831 CHKERRQ(MatMissingDiagonal(B,&miss,&d)); 832 CHKERRQ(MatGetOwnershipRange(B,&rst,&ren)); 833 CHKERRQ(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 834 CHKERRQ(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "[%d] [%" PetscInt_FMT ",%" PetscInt_FMT ") Missing %d, row %" PetscInt_FMT " (diag %s)\n",rank,rst,ren,(int)miss,d,diagstr)); 835 CHKERRQ(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 836 CHKERRQ(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD)); 837 838 CHKERRQ(VecDestroy(&x)); 839 CHKERRQ(VecDestroy(&b)); 840 CHKERRQ(VecDestroy(&b2)); 841 842 /* check the result of ZeroRows with that from MPIAIJ routines 843 assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */ 844 if (haszerorows) { 845 CHKERRQ(MatDuplicate(Afull,MAT_COPY_VALUES,&Bcheck)); 846 CHKERRQ(MatSetOption(Bcheck,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 847 CHKERRQ(MatZeroRowsIS(Bcheck,is,diag,NULL,NULL)); 848 CHKERRQ(CheckMat(B,Bcheck,PETSC_FALSE,"Zerorows")); 849 CHKERRQ(MatDestroy(&Bcheck)); 850 } 851 CHKERRQ(MatDestroy(&B)); 852 853 if (B2) { /* test MatZeroRowsColumns */ 854 CHKERRQ(MatDuplicate(Afull,MAT_COPY_VALUES,&B)); 855 CHKERRQ(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 856 CHKERRQ(MatZeroRowsColumnsIS(B,is,diag,NULL,NULL)); 857 CHKERRQ(CheckMat(B2,B,PETSC_FALSE,"MatZeroRowsColumns")); 858 CHKERRQ(MatDestroy(&B)); 859 CHKERRQ(MatDestroy(&B2)); 860 } 861 PetscFunctionReturn(0); 862 } 863 864 /*TEST 865 866 test: 867 args: -test_trans 868 869 test: 870 suffix: 2 871 nsize: 4 872 args: -matis_convert_local_nest -nr 3 -nc 4 873 874 test: 875 suffix: 3 876 nsize: 5 877 args: -m 11 -n 10 -matis_convert_local_nest -nr 2 -nc 1 878 879 test: 880 suffix: 4 881 nsize: 6 882 args: -m 9 -n 12 -test_trans -nr 2 -nc 7 883 884 test: 885 suffix: 5 886 nsize: 6 887 args: -m 12 -n 12 -test_trans -nr 3 -nc 1 888 889 test: 890 suffix: 6 891 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap 892 893 test: 894 suffix: 7 895 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap 896 897 test: 898 suffix: 8 899 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap 900 901 test: 902 suffix: 9 903 nsize: 5 904 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap 905 906 test: 907 suffix: 10 908 nsize: 5 909 args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap 910 911 test: 912 suffix: vscat_default 913 nsize: 5 914 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap 915 output_file: output/ex23_11.out 916 917 test: 918 suffix: 12 919 nsize: 3 920 args: -m 12 -n 12 -symmetric -matis_localmat_type sbaij -test_trans -nr 2 -nc 3 921 922 testset: 923 output_file: output/ex23_13.out 924 nsize: 3 925 args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap 926 filter: grep -v "type:" 927 test: 928 suffix: baij 929 args: -matis_localmat_type baij 930 test: 931 requires: viennacl 932 suffix: viennacl 933 args: -matis_localmat_type aijviennacl 934 test: 935 requires: cuda 936 suffix: cusparse 937 args: -matis_localmat_type aijcusparse 938 939 test: 940 suffix: negrep 941 nsize: {{1 3}separate output} 942 args: -m {{5 7}separate output} -n {{5 7}separate output} -test_trans -nr 2 -nc 3 -negmap {{0 1}separate output} -repmap {{0 1}separate output} -permmap -diffmap {{0 1}separate output} 943 944 TEST*/ 945