1 2 #include <../src/mat/impls/baij/mpi/mpibaij.h> /*I "petscmat.h" I*/ 3 4 #include <petscblaslapack.h> 5 #include <petscsf.h> 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "MatGetRowMaxAbs_MPIBAIJ" 9 PetscErrorCode MatGetRowMaxAbs_MPIBAIJ(Mat A,Vec v,PetscInt idx[]) 10 { 11 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 12 PetscErrorCode ierr; 13 PetscInt i,*idxb = 0; 14 PetscScalar *va,*vb; 15 Vec vtmp; 16 17 PetscFunctionBegin; 18 ierr = MatGetRowMaxAbs(a->A,v,idx);CHKERRQ(ierr); 19 ierr = VecGetArray(v,&va);CHKERRQ(ierr); 20 if (idx) { 21 for (i=0; i<A->rmap->n; i++) { 22 if (PetscAbsScalar(va[i])) idx[i] += A->cmap->rstart; 23 } 24 } 25 26 ierr = VecCreateSeq(PETSC_COMM_SELF,A->rmap->n,&vtmp);CHKERRQ(ierr); 27 if (idx) {ierr = PetscMalloc1(A->rmap->n,&idxb);CHKERRQ(ierr);} 28 ierr = MatGetRowMaxAbs(a->B,vtmp,idxb);CHKERRQ(ierr); 29 ierr = VecGetArray(vtmp,&vb);CHKERRQ(ierr); 30 31 for (i=0; i<A->rmap->n; i++) { 32 if (PetscAbsScalar(va[i]) < PetscAbsScalar(vb[i])) { 33 va[i] = vb[i]; 34 if (idx) idx[i] = A->cmap->bs*a->garray[idxb[i]/A->cmap->bs] + (idxb[i] % A->cmap->bs); 35 } 36 } 37 38 ierr = VecRestoreArray(v,&va);CHKERRQ(ierr); 39 ierr = VecRestoreArray(vtmp,&vb);CHKERRQ(ierr); 40 ierr = PetscFree(idxb);CHKERRQ(ierr); 41 ierr = VecDestroy(&vtmp);CHKERRQ(ierr); 42 PetscFunctionReturn(0); 43 } 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "MatStoreValues_MPIBAIJ" 47 PetscErrorCode MatStoreValues_MPIBAIJ(Mat mat) 48 { 49 Mat_MPIBAIJ *aij = (Mat_MPIBAIJ*)mat->data; 50 PetscErrorCode ierr; 51 52 PetscFunctionBegin; 53 ierr = MatStoreValues(aij->A);CHKERRQ(ierr); 54 ierr = MatStoreValues(aij->B);CHKERRQ(ierr); 55 PetscFunctionReturn(0); 56 } 57 58 #undef __FUNCT__ 59 #define __FUNCT__ "MatRetrieveValues_MPIBAIJ" 60 PetscErrorCode MatRetrieveValues_MPIBAIJ(Mat mat) 61 { 62 Mat_MPIBAIJ *aij = (Mat_MPIBAIJ*)mat->data; 63 PetscErrorCode ierr; 64 65 PetscFunctionBegin; 66 ierr = MatRetrieveValues(aij->A);CHKERRQ(ierr); 67 ierr = MatRetrieveValues(aij->B);CHKERRQ(ierr); 68 PetscFunctionReturn(0); 69 } 70 71 /* 72 Local utility routine that creates a mapping from the global column 73 number to the local number in the off-diagonal part of the local 74 storage of the matrix. This is done in a non scalable way since the 75 length of colmap equals the global matrix length. 76 */ 77 #undef __FUNCT__ 78 #define __FUNCT__ "MatCreateColmap_MPIBAIJ_Private" 79 PetscErrorCode MatCreateColmap_MPIBAIJ_Private(Mat mat) 80 { 81 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 82 Mat_SeqBAIJ *B = (Mat_SeqBAIJ*)baij->B->data; 83 PetscErrorCode ierr; 84 PetscInt nbs = B->nbs,i,bs=mat->rmap->bs; 85 86 PetscFunctionBegin; 87 #if defined(PETSC_USE_CTABLE) 88 ierr = PetscTableCreate(baij->nbs,baij->Nbs+1,&baij->colmap);CHKERRQ(ierr); 89 for (i=0; i<nbs; i++) { 90 ierr = PetscTableAdd(baij->colmap,baij->garray[i]+1,i*bs+1,INSERT_VALUES);CHKERRQ(ierr); 91 } 92 #else 93 ierr = PetscMalloc1(baij->Nbs+1,&baij->colmap);CHKERRQ(ierr); 94 ierr = PetscLogObjectMemory((PetscObject)mat,baij->Nbs*sizeof(PetscInt));CHKERRQ(ierr); 95 ierr = PetscMemzero(baij->colmap,baij->Nbs*sizeof(PetscInt));CHKERRQ(ierr); 96 for (i=0; i<nbs; i++) baij->colmap[baij->garray[i]] = i*bs+1; 97 #endif 98 PetscFunctionReturn(0); 99 } 100 101 #define MatSetValues_SeqBAIJ_A_Private(row,col,value,addv,orow,ocol) \ 102 { \ 103 \ 104 brow = row/bs; \ 105 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; \ 106 rmax = aimax[brow]; nrow = ailen[brow]; \ 107 bcol = col/bs; \ 108 ridx = row % bs; cidx = col % bs; \ 109 low = 0; high = nrow; \ 110 while (high-low > 3) { \ 111 t = (low+high)/2; \ 112 if (rp[t] > bcol) high = t; \ 113 else low = t; \ 114 } \ 115 for (_i=low; _i<high; _i++) { \ 116 if (rp[_i] > bcol) break; \ 117 if (rp[_i] == bcol) { \ 118 bap = ap + bs2*_i + bs*cidx + ridx; \ 119 if (addv == ADD_VALUES) *bap += value; \ 120 else *bap = value; \ 121 goto a_noinsert; \ 122 } \ 123 } \ 124 if (a->nonew == 1) goto a_noinsert; \ 125 if (a->nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%D, %D) into matrix", orow, ocol); \ 126 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,aimax,a->nonew,MatScalar); \ 127 N = nrow++ - 1; \ 128 /* shift up all the later entries in this row */ \ 129 for (ii=N; ii>=_i; ii--) { \ 130 rp[ii+1] = rp[ii]; \ 131 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); \ 132 } \ 133 if (N>=_i) { ierr = PetscMemzero(ap+bs2*_i,bs2*sizeof(MatScalar));CHKERRQ(ierr); } \ 134 rp[_i] = bcol; \ 135 ap[bs2*_i + bs*cidx + ridx] = value; \ 136 a_noinsert:; \ 137 ailen[brow] = nrow; \ 138 } 139 140 #define MatSetValues_SeqBAIJ_B_Private(row,col,value,addv,orow,ocol) \ 141 { \ 142 brow = row/bs; \ 143 rp = bj + bi[brow]; ap = ba + bs2*bi[brow]; \ 144 rmax = bimax[brow]; nrow = bilen[brow]; \ 145 bcol = col/bs; \ 146 ridx = row % bs; cidx = col % bs; \ 147 low = 0; high = nrow; \ 148 while (high-low > 3) { \ 149 t = (low+high)/2; \ 150 if (rp[t] > bcol) high = t; \ 151 else low = t; \ 152 } \ 153 for (_i=low; _i<high; _i++) { \ 154 if (rp[_i] > bcol) break; \ 155 if (rp[_i] == bcol) { \ 156 bap = ap + bs2*_i + bs*cidx + ridx; \ 157 if (addv == ADD_VALUES) *bap += value; \ 158 else *bap = value; \ 159 goto b_noinsert; \ 160 } \ 161 } \ 162 if (b->nonew == 1) goto b_noinsert; \ 163 if (b->nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%D, %D) into matrix", orow, ocol); \ 164 MatSeqXAIJReallocateAIJ(B,b->mbs,bs2,nrow,brow,bcol,rmax,ba,bi,bj,rp,ap,bimax,b->nonew,MatScalar); \ 165 N = nrow++ - 1; \ 166 /* shift up all the later entries in this row */ \ 167 for (ii=N; ii>=_i; ii--) { \ 168 rp[ii+1] = rp[ii]; \ 169 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); \ 170 } \ 171 if (N>=_i) { ierr = PetscMemzero(ap+bs2*_i,bs2*sizeof(MatScalar));CHKERRQ(ierr);} \ 172 rp[_i] = bcol; \ 173 ap[bs2*_i + bs*cidx + ridx] = value; \ 174 b_noinsert:; \ 175 bilen[brow] = nrow; \ 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "MatSetValues_MPIBAIJ" 180 PetscErrorCode MatSetValues_MPIBAIJ(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv) 181 { 182 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 183 MatScalar value; 184 PetscBool roworiented = baij->roworiented; 185 PetscErrorCode ierr; 186 PetscInt i,j,row,col; 187 PetscInt rstart_orig=mat->rmap->rstart; 188 PetscInt rend_orig =mat->rmap->rend,cstart_orig=mat->cmap->rstart; 189 PetscInt cend_orig =mat->cmap->rend,bs=mat->rmap->bs; 190 191 /* Some Variables required in the macro */ 192 Mat A = baij->A; 193 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)(A)->data; 194 PetscInt *aimax=a->imax,*ai=a->i,*ailen=a->ilen,*aj=a->j; 195 MatScalar *aa =a->a; 196 197 Mat B = baij->B; 198 Mat_SeqBAIJ *b = (Mat_SeqBAIJ*)(B)->data; 199 PetscInt *bimax=b->imax,*bi=b->i,*bilen=b->ilen,*bj=b->j; 200 MatScalar *ba =b->a; 201 202 PetscInt *rp,ii,nrow,_i,rmax,N,brow,bcol; 203 PetscInt low,high,t,ridx,cidx,bs2=a->bs2; 204 MatScalar *ap,*bap; 205 206 PetscFunctionBegin; 207 for (i=0; i<m; i++) { 208 if (im[i] < 0) continue; 209 #if defined(PETSC_USE_DEBUG) 210 if (im[i] >= mat->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",im[i],mat->rmap->N-1); 211 #endif 212 if (im[i] >= rstart_orig && im[i] < rend_orig) { 213 row = im[i] - rstart_orig; 214 for (j=0; j<n; j++) { 215 if (in[j] >= cstart_orig && in[j] < cend_orig) { 216 col = in[j] - cstart_orig; 217 if (roworiented) value = v[i*n+j]; 218 else value = v[i+j*m]; 219 MatSetValues_SeqBAIJ_A_Private(row,col,value,addv,im[i],in[j]); 220 /* ierr = MatSetValues_SeqBAIJ(baij->A,1,&row,1,&col,&value,addv);CHKERRQ(ierr); */ 221 } else if (in[j] < 0) continue; 222 #if defined(PETSC_USE_DEBUG) 223 else if (in[j] >= mat->cmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[j],mat->cmap->N-1); 224 #endif 225 else { 226 if (mat->was_assembled) { 227 if (!baij->colmap) { 228 ierr = MatCreateColmap_MPIBAIJ_Private(mat);CHKERRQ(ierr); 229 } 230 #if defined(PETSC_USE_CTABLE) 231 ierr = PetscTableFind(baij->colmap,in[j]/bs + 1,&col);CHKERRQ(ierr); 232 col = col - 1; 233 #else 234 col = baij->colmap[in[j]/bs] - 1; 235 #endif 236 if (col < 0 && !((Mat_SeqBAIJ*)(baij->B->data))->nonew) { 237 ierr = MatDisAssemble_MPIBAIJ(mat);CHKERRQ(ierr); 238 col = in[j]; 239 /* Reinitialize the variables required by MatSetValues_SeqBAIJ_B_Private() */ 240 B = baij->B; 241 b = (Mat_SeqBAIJ*)(B)->data; 242 bimax=b->imax;bi=b->i;bilen=b->ilen;bj=b->j; 243 ba =b->a; 244 } else if (col < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) into matrix", im[i], in[j]); 245 else col += in[j]%bs; 246 } else col = in[j]; 247 if (roworiented) value = v[i*n+j]; 248 else value = v[i+j*m]; 249 MatSetValues_SeqBAIJ_B_Private(row,col,value,addv,im[i],in[j]); 250 /* ierr = MatSetValues_SeqBAIJ(baij->B,1,&row,1,&col,&value,addv);CHKERRQ(ierr); */ 251 } 252 } 253 } else { 254 if (mat->nooffprocentries) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Setting off process row %D even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set",im[i]); 255 if (!baij->donotstash) { 256 mat->assembled = PETSC_FALSE; 257 if (roworiented) { 258 ierr = MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n,PETSC_FALSE);CHKERRQ(ierr); 259 } else { 260 ierr = MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m,PETSC_FALSE);CHKERRQ(ierr); 261 } 262 } 263 } 264 } 265 PetscFunctionReturn(0); 266 } 267 268 #undef __FUNCT__ 269 #define __FUNCT__ "MatSetValuesBlocked_SeqBAIJ_Inlined" 270 PETSC_STATIC_INLINE PetscErrorCode MatSetValuesBlocked_SeqBAIJ_Inlined(Mat A,PetscInt row,PetscInt col,const PetscScalar v[],InsertMode is,PetscInt orow,PetscInt ocol) 271 { 272 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 273 PetscInt *rp,low,high,t,ii,jj,nrow,i,rmax,N; 274 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen; 275 PetscErrorCode ierr; 276 PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs; 277 PetscBool roworiented=a->roworiented; 278 const PetscScalar *value = v; 279 MatScalar *ap,*aa = a->a,*bap; 280 281 PetscFunctionBegin; 282 rp = aj + ai[row]; 283 ap = aa + bs2*ai[row]; 284 rmax = imax[row]; 285 nrow = ailen[row]; 286 value = v; 287 low = 0; 288 high = nrow; 289 while (high-low > 7) { 290 t = (low+high)/2; 291 if (rp[t] > col) high = t; 292 else low = t; 293 } 294 for (i=low; i<high; i++) { 295 if (rp[i] > col) break; 296 if (rp[i] == col) { 297 bap = ap + bs2*i; 298 if (roworiented) { 299 if (is == ADD_VALUES) { 300 for (ii=0; ii<bs; ii++) { 301 for (jj=ii; jj<bs2; jj+=bs) { 302 bap[jj] += *value++; 303 } 304 } 305 } else { 306 for (ii=0; ii<bs; ii++) { 307 for (jj=ii; jj<bs2; jj+=bs) { 308 bap[jj] = *value++; 309 } 310 } 311 } 312 } else { 313 if (is == ADD_VALUES) { 314 for (ii=0; ii<bs; ii++,value+=bs) { 315 for (jj=0; jj<bs; jj++) { 316 bap[jj] += value[jj]; 317 } 318 bap += bs; 319 } 320 } else { 321 for (ii=0; ii<bs; ii++,value+=bs) { 322 for (jj=0; jj<bs; jj++) { 323 bap[jj] = value[jj]; 324 } 325 bap += bs; 326 } 327 } 328 } 329 goto noinsert2; 330 } 331 } 332 if (nonew == 1) goto noinsert2; 333 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new global block indexed nonzero block (%D, %D) in the matrix", orow, ocol); 334 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 335 N = nrow++ - 1; high++; 336 /* shift up all the later entries in this row */ 337 for (ii=N; ii>=i; ii--) { 338 rp[ii+1] = rp[ii]; 339 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 340 } 341 if (N >= i) { 342 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 343 } 344 rp[i] = col; 345 bap = ap + bs2*i; 346 if (roworiented) { 347 for (ii=0; ii<bs; ii++) { 348 for (jj=ii; jj<bs2; jj+=bs) { 349 bap[jj] = *value++; 350 } 351 } 352 } else { 353 for (ii=0; ii<bs; ii++) { 354 for (jj=0; jj<bs; jj++) { 355 *bap++ = *value++; 356 } 357 } 358 } 359 noinsert2:; 360 ailen[row] = nrow; 361 PetscFunctionReturn(0); 362 } 363 364 #undef __FUNCT__ 365 #define __FUNCT__ "MatSetValuesBlocked_MPIBAIJ" 366 /* 367 This routine should be optimized so that the block copy at ** Here a copy is required ** below is not needed 368 by passing additional stride information into the MatSetValuesBlocked_SeqBAIJ_Inlined() routine 369 */ 370 PetscErrorCode MatSetValuesBlocked_MPIBAIJ(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv) 371 { 372 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 373 const PetscScalar *value; 374 MatScalar *barray = baij->barray; 375 PetscBool roworiented = baij->roworiented; 376 PetscErrorCode ierr; 377 PetscInt i,j,ii,jj,row,col,rstart=baij->rstartbs; 378 PetscInt rend=baij->rendbs,cstart=baij->cstartbs,stepval; 379 PetscInt cend=baij->cendbs,bs=mat->rmap->bs,bs2=baij->bs2; 380 381 PetscFunctionBegin; 382 if (!barray) { 383 ierr = PetscMalloc1(bs2,&barray);CHKERRQ(ierr); 384 baij->barray = barray; 385 } 386 387 if (roworiented) stepval = (n-1)*bs; 388 else stepval = (m-1)*bs; 389 390 for (i=0; i<m; i++) { 391 if (im[i] < 0) continue; 392 #if defined(PETSC_USE_DEBUG) 393 if (im[i] >= baij->Mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block indexed row too large %D max %D",im[i],baij->Mbs-1); 394 #endif 395 if (im[i] >= rstart && im[i] < rend) { 396 row = im[i] - rstart; 397 for (j=0; j<n; j++) { 398 /* If NumCol = 1 then a copy is not required */ 399 if ((roworiented) && (n == 1)) { 400 barray = (MatScalar*)v + i*bs2; 401 } else if ((!roworiented) && (m == 1)) { 402 barray = (MatScalar*)v + j*bs2; 403 } else { /* Here a copy is required */ 404 if (roworiented) { 405 value = v + (i*(stepval+bs) + j)*bs; 406 } else { 407 value = v + (j*(stepval+bs) + i)*bs; 408 } 409 for (ii=0; ii<bs; ii++,value+=bs+stepval) { 410 for (jj=0; jj<bs; jj++) barray[jj] = value[jj]; 411 barray += bs; 412 } 413 barray -= bs2; 414 } 415 416 if (in[j] >= cstart && in[j] < cend) { 417 col = in[j] - cstart; 418 ierr = MatSetValuesBlocked_SeqBAIJ_Inlined(baij->A,row,col,barray,addv,im[i],in[j]);CHKERRQ(ierr); 419 } else if (in[j] < 0) continue; 420 #if defined(PETSC_USE_DEBUG) 421 else if (in[j] >= baij->Nbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block indexed column too large %D max %D",in[j],baij->Nbs-1); 422 #endif 423 else { 424 if (mat->was_assembled) { 425 if (!baij->colmap) { 426 ierr = MatCreateColmap_MPIBAIJ_Private(mat);CHKERRQ(ierr); 427 } 428 429 #if defined(PETSC_USE_DEBUG) 430 #if defined(PETSC_USE_CTABLE) 431 { PetscInt data; 432 ierr = PetscTableFind(baij->colmap,in[j]+1,&data);CHKERRQ(ierr); 433 if ((data - 1) % bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Incorrect colmap"); 434 } 435 #else 436 if ((baij->colmap[in[j]] - 1) % bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Incorrect colmap"); 437 #endif 438 #endif 439 #if defined(PETSC_USE_CTABLE) 440 ierr = PetscTableFind(baij->colmap,in[j]+1,&col);CHKERRQ(ierr); 441 col = (col - 1)/bs; 442 #else 443 col = (baij->colmap[in[j]] - 1)/bs; 444 #endif 445 if (col < 0 && !((Mat_SeqBAIJ*)(baij->B->data))->nonew) { 446 ierr = MatDisAssemble_MPIBAIJ(mat);CHKERRQ(ierr); 447 col = in[j]; 448 } else if (col < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new blocked indexed nonzero block (%D, %D) into matrix",im[i],in[j]); 449 } else col = in[j]; 450 ierr = MatSetValuesBlocked_SeqBAIJ_Inlined(baij->B,row,col,barray,addv,im[i],in[j]);CHKERRQ(ierr); 451 } 452 } 453 } else { 454 if (mat->nooffprocentries) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Setting off process block indexed row %D even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set",im[i]); 455 if (!baij->donotstash) { 456 if (roworiented) { 457 ierr = MatStashValuesRowBlocked_Private(&mat->bstash,im[i],n,in,v,m,n,i);CHKERRQ(ierr); 458 } else { 459 ierr = MatStashValuesColBlocked_Private(&mat->bstash,im[i],n,in,v,m,n,i);CHKERRQ(ierr); 460 } 461 } 462 } 463 } 464 PetscFunctionReturn(0); 465 } 466 467 #define HASH_KEY 0.6180339887 468 #define HASH(size,key,tmp) (tmp = (key)*HASH_KEY,(PetscInt)((size)*(tmp-(PetscInt)tmp))) 469 /* #define HASH(size,key) ((PetscInt)((size)*fmod(((key)*HASH_KEY),1))) */ 470 /* #define HASH(size,key,tmp) ((PetscInt)((size)*fmod(((key)*HASH_KEY),1))) */ 471 #undef __FUNCT__ 472 #define __FUNCT__ "MatSetValues_MPIBAIJ_HT" 473 PetscErrorCode MatSetValues_MPIBAIJ_HT(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv) 474 { 475 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 476 PetscBool roworiented = baij->roworiented; 477 PetscErrorCode ierr; 478 PetscInt i,j,row,col; 479 PetscInt rstart_orig=mat->rmap->rstart; 480 PetscInt rend_orig =mat->rmap->rend,Nbs=baij->Nbs; 481 PetscInt h1,key,size=baij->ht_size,bs=mat->rmap->bs,*HT=baij->ht,idx; 482 PetscReal tmp; 483 MatScalar **HD = baij->hd,value; 484 #if defined(PETSC_USE_DEBUG) 485 PetscInt total_ct=baij->ht_total_ct,insert_ct=baij->ht_insert_ct; 486 #endif 487 488 PetscFunctionBegin; 489 for (i=0; i<m; i++) { 490 #if defined(PETSC_USE_DEBUG) 491 if (im[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row"); 492 if (im[i] >= mat->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",im[i],mat->rmap->N-1); 493 #endif 494 row = im[i]; 495 if (row >= rstart_orig && row < rend_orig) { 496 for (j=0; j<n; j++) { 497 col = in[j]; 498 if (roworiented) value = v[i*n+j]; 499 else value = v[i+j*m]; 500 /* Look up PetscInto the Hash Table */ 501 key = (row/bs)*Nbs+(col/bs)+1; 502 h1 = HASH(size,key,tmp); 503 504 505 idx = h1; 506 #if defined(PETSC_USE_DEBUG) 507 insert_ct++; 508 total_ct++; 509 if (HT[idx] != key) { 510 for (idx=h1; (idx<size) && (HT[idx]!=key); idx++,total_ct++) ; 511 if (idx == size) { 512 for (idx=0; (idx<h1) && (HT[idx]!=key); idx++,total_ct++) ; 513 if (idx == h1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"(%D,%D) has no entry in the hash table", row, col); 514 } 515 } 516 #else 517 if (HT[idx] != key) { 518 for (idx=h1; (idx<size) && (HT[idx]!=key); idx++) ; 519 if (idx == size) { 520 for (idx=0; (idx<h1) && (HT[idx]!=key); idx++) ; 521 if (idx == h1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"(%D,%D) has no entry in the hash table", row, col); 522 } 523 } 524 #endif 525 /* A HASH table entry is found, so insert the values at the correct address */ 526 if (addv == ADD_VALUES) *(HD[idx]+ (col % bs)*bs + (row % bs)) += value; 527 else *(HD[idx]+ (col % bs)*bs + (row % bs)) = value; 528 } 529 } else if (!baij->donotstash) { 530 if (roworiented) { 531 ierr = MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n,PETSC_FALSE);CHKERRQ(ierr); 532 } else { 533 ierr = MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m,PETSC_FALSE);CHKERRQ(ierr); 534 } 535 } 536 } 537 #if defined(PETSC_USE_DEBUG) 538 baij->ht_total_ct += total_ct; 539 baij->ht_insert_ct += insert_ct; 540 #endif 541 PetscFunctionReturn(0); 542 } 543 544 #undef __FUNCT__ 545 #define __FUNCT__ "MatSetValuesBlocked_MPIBAIJ_HT" 546 PetscErrorCode MatSetValuesBlocked_MPIBAIJ_HT(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv) 547 { 548 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 549 PetscBool roworiented = baij->roworiented; 550 PetscErrorCode ierr; 551 PetscInt i,j,ii,jj,row,col; 552 PetscInt rstart=baij->rstartbs; 553 PetscInt rend =mat->rmap->rend,stepval,bs=mat->rmap->bs,bs2=baij->bs2,nbs2=n*bs2; 554 PetscInt h1,key,size=baij->ht_size,idx,*HT=baij->ht,Nbs=baij->Nbs; 555 PetscReal tmp; 556 MatScalar **HD = baij->hd,*baij_a; 557 const PetscScalar *v_t,*value; 558 #if defined(PETSC_USE_DEBUG) 559 PetscInt total_ct=baij->ht_total_ct,insert_ct=baij->ht_insert_ct; 560 #endif 561 562 PetscFunctionBegin; 563 if (roworiented) stepval = (n-1)*bs; 564 else stepval = (m-1)*bs; 565 566 for (i=0; i<m; i++) { 567 #if defined(PETSC_USE_DEBUG) 568 if (im[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",im[i]); 569 if (im[i] >= baij->Mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",im[i],baij->Mbs-1); 570 #endif 571 row = im[i]; 572 v_t = v + i*nbs2; 573 if (row >= rstart && row < rend) { 574 for (j=0; j<n; j++) { 575 col = in[j]; 576 577 /* Look up into the Hash Table */ 578 key = row*Nbs+col+1; 579 h1 = HASH(size,key,tmp); 580 581 idx = h1; 582 #if defined(PETSC_USE_DEBUG) 583 total_ct++; 584 insert_ct++; 585 if (HT[idx] != key) { 586 for (idx=h1; (idx<size) && (HT[idx]!=key); idx++,total_ct++) ; 587 if (idx == size) { 588 for (idx=0; (idx<h1) && (HT[idx]!=key); idx++,total_ct++) ; 589 if (idx == h1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"(%D,%D) has no entry in the hash table", row, col); 590 } 591 } 592 #else 593 if (HT[idx] != key) { 594 for (idx=h1; (idx<size) && (HT[idx]!=key); idx++) ; 595 if (idx == size) { 596 for (idx=0; (idx<h1) && (HT[idx]!=key); idx++) ; 597 if (idx == h1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"(%D,%D) has no entry in the hash table", row, col); 598 } 599 } 600 #endif 601 baij_a = HD[idx]; 602 if (roworiented) { 603 /*value = v + i*(stepval+bs)*bs + j*bs;*/ 604 /* value = v + (i*(stepval+bs)+j)*bs; */ 605 value = v_t; 606 v_t += bs; 607 if (addv == ADD_VALUES) { 608 for (ii=0; ii<bs; ii++,value+=stepval) { 609 for (jj=ii; jj<bs2; jj+=bs) { 610 baij_a[jj] += *value++; 611 } 612 } 613 } else { 614 for (ii=0; ii<bs; ii++,value+=stepval) { 615 for (jj=ii; jj<bs2; jj+=bs) { 616 baij_a[jj] = *value++; 617 } 618 } 619 } 620 } else { 621 value = v + j*(stepval+bs)*bs + i*bs; 622 if (addv == ADD_VALUES) { 623 for (ii=0; ii<bs; ii++,value+=stepval,baij_a+=bs) { 624 for (jj=0; jj<bs; jj++) { 625 baij_a[jj] += *value++; 626 } 627 } 628 } else { 629 for (ii=0; ii<bs; ii++,value+=stepval,baij_a+=bs) { 630 for (jj=0; jj<bs; jj++) { 631 baij_a[jj] = *value++; 632 } 633 } 634 } 635 } 636 } 637 } else { 638 if (!baij->donotstash) { 639 if (roworiented) { 640 ierr = MatStashValuesRowBlocked_Private(&mat->bstash,im[i],n,in,v,m,n,i);CHKERRQ(ierr); 641 } else { 642 ierr = MatStashValuesColBlocked_Private(&mat->bstash,im[i],n,in,v,m,n,i);CHKERRQ(ierr); 643 } 644 } 645 } 646 } 647 #if defined(PETSC_USE_DEBUG) 648 baij->ht_total_ct += total_ct; 649 baij->ht_insert_ct += insert_ct; 650 #endif 651 PetscFunctionReturn(0); 652 } 653 654 #undef __FUNCT__ 655 #define __FUNCT__ "MatGetValues_MPIBAIJ" 656 PetscErrorCode MatGetValues_MPIBAIJ(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[]) 657 { 658 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 659 PetscErrorCode ierr; 660 PetscInt bs = mat->rmap->bs,i,j,bsrstart = mat->rmap->rstart,bsrend = mat->rmap->rend; 661 PetscInt bscstart = mat->cmap->rstart,bscend = mat->cmap->rend,row,col,data; 662 663 PetscFunctionBegin; 664 for (i=0; i<m; i++) { 665 if (idxm[i] < 0) continue; /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",idxm[i]);*/ 666 if (idxm[i] >= mat->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm[i],mat->rmap->N-1); 667 if (idxm[i] >= bsrstart && idxm[i] < bsrend) { 668 row = idxm[i] - bsrstart; 669 for (j=0; j<n; j++) { 670 if (idxn[j] < 0) continue; /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",idxn[j]); */ 671 if (idxn[j] >= mat->cmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",idxn[j],mat->cmap->N-1); 672 if (idxn[j] >= bscstart && idxn[j] < bscend) { 673 col = idxn[j] - bscstart; 674 ierr = MatGetValues_SeqBAIJ(baij->A,1,&row,1,&col,v+i*n+j);CHKERRQ(ierr); 675 } else { 676 if (!baij->colmap) { 677 ierr = MatCreateColmap_MPIBAIJ_Private(mat);CHKERRQ(ierr); 678 } 679 #if defined(PETSC_USE_CTABLE) 680 ierr = PetscTableFind(baij->colmap,idxn[j]/bs+1,&data);CHKERRQ(ierr); 681 data--; 682 #else 683 data = baij->colmap[idxn[j]/bs]-1; 684 #endif 685 if ((data < 0) || (baij->garray[data/bs] != idxn[j]/bs)) *(v+i*n+j) = 0.0; 686 else { 687 col = data + idxn[j]%bs; 688 ierr = MatGetValues_SeqBAIJ(baij->B,1,&row,1,&col,v+i*n+j);CHKERRQ(ierr); 689 } 690 } 691 } 692 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only local values currently supported"); 693 } 694 PetscFunctionReturn(0); 695 } 696 697 #undef __FUNCT__ 698 #define __FUNCT__ "MatNorm_MPIBAIJ" 699 PetscErrorCode MatNorm_MPIBAIJ(Mat mat,NormType type,PetscReal *nrm) 700 { 701 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 702 Mat_SeqBAIJ *amat = (Mat_SeqBAIJ*)baij->A->data,*bmat = (Mat_SeqBAIJ*)baij->B->data; 703 PetscErrorCode ierr; 704 PetscInt i,j,bs2=baij->bs2,bs=baij->A->rmap->bs,nz,row,col; 705 PetscReal sum = 0.0; 706 MatScalar *v; 707 708 PetscFunctionBegin; 709 if (baij->size == 1) { 710 ierr = MatNorm(baij->A,type,nrm);CHKERRQ(ierr); 711 } else { 712 if (type == NORM_FROBENIUS) { 713 v = amat->a; 714 nz = amat->nz*bs2; 715 for (i=0; i<nz; i++) { 716 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 717 } 718 v = bmat->a; 719 nz = bmat->nz*bs2; 720 for (i=0; i<nz; i++) { 721 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 722 } 723 ierr = MPIU_Allreduce(&sum,nrm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 724 *nrm = PetscSqrtReal(*nrm); 725 } else if (type == NORM_1) { /* max column sum */ 726 PetscReal *tmp,*tmp2; 727 PetscInt *jj,*garray=baij->garray,cstart=baij->rstartbs; 728 ierr = PetscMalloc2(mat->cmap->N,&tmp,mat->cmap->N,&tmp2);CHKERRQ(ierr); 729 ierr = PetscMemzero(tmp,mat->cmap->N*sizeof(PetscReal));CHKERRQ(ierr); 730 v = amat->a; jj = amat->j; 731 for (i=0; i<amat->nz; i++) { 732 for (j=0; j<bs; j++) { 733 col = bs*(cstart + *jj) + j; /* column index */ 734 for (row=0; row<bs; row++) { 735 tmp[col] += PetscAbsScalar(*v); v++; 736 } 737 } 738 jj++; 739 } 740 v = bmat->a; jj = bmat->j; 741 for (i=0; i<bmat->nz; i++) { 742 for (j=0; j<bs; j++) { 743 col = bs*garray[*jj] + j; 744 for (row=0; row<bs; row++) { 745 tmp[col] += PetscAbsScalar(*v); v++; 746 } 747 } 748 jj++; 749 } 750 ierr = MPIU_Allreduce(tmp,tmp2,mat->cmap->N,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 751 *nrm = 0.0; 752 for (j=0; j<mat->cmap->N; j++) { 753 if (tmp2[j] > *nrm) *nrm = tmp2[j]; 754 } 755 ierr = PetscFree2(tmp,tmp2);CHKERRQ(ierr); 756 } else if (type == NORM_INFINITY) { /* max row sum */ 757 PetscReal *sums; 758 ierr = PetscMalloc1(bs,&sums);CHKERRQ(ierr); 759 sum = 0.0; 760 for (j=0; j<amat->mbs; j++) { 761 for (row=0; row<bs; row++) sums[row] = 0.0; 762 v = amat->a + bs2*amat->i[j]; 763 nz = amat->i[j+1]-amat->i[j]; 764 for (i=0; i<nz; i++) { 765 for (col=0; col<bs; col++) { 766 for (row=0; row<bs; row++) { 767 sums[row] += PetscAbsScalar(*v); v++; 768 } 769 } 770 } 771 v = bmat->a + bs2*bmat->i[j]; 772 nz = bmat->i[j+1]-bmat->i[j]; 773 for (i=0; i<nz; i++) { 774 for (col=0; col<bs; col++) { 775 for (row=0; row<bs; row++) { 776 sums[row] += PetscAbsScalar(*v); v++; 777 } 778 } 779 } 780 for (row=0; row<bs; row++) { 781 if (sums[row] > sum) sum = sums[row]; 782 } 783 } 784 ierr = MPIU_Allreduce(&sum,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 785 ierr = PetscFree(sums);CHKERRQ(ierr); 786 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support for this norm yet"); 787 } 788 PetscFunctionReturn(0); 789 } 790 791 /* 792 Creates the hash table, and sets the table 793 This table is created only once. 794 If new entried need to be added to the matrix 795 then the hash table has to be destroyed and 796 recreated. 797 */ 798 #undef __FUNCT__ 799 #define __FUNCT__ "MatCreateHashTable_MPIBAIJ_Private" 800 PetscErrorCode MatCreateHashTable_MPIBAIJ_Private(Mat mat,PetscReal factor) 801 { 802 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 803 Mat A = baij->A,B=baij->B; 804 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)B->data; 805 PetscInt i,j,k,nz=a->nz+b->nz,h1,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j; 806 PetscErrorCode ierr; 807 PetscInt ht_size,bs2=baij->bs2,rstart=baij->rstartbs; 808 PetscInt cstart=baij->cstartbs,*garray=baij->garray,row,col,Nbs=baij->Nbs; 809 PetscInt *HT,key; 810 MatScalar **HD; 811 PetscReal tmp; 812 #if defined(PETSC_USE_INFO) 813 PetscInt ct=0,max=0; 814 #endif 815 816 PetscFunctionBegin; 817 if (baij->ht) PetscFunctionReturn(0); 818 819 baij->ht_size = (PetscInt)(factor*nz); 820 ht_size = baij->ht_size; 821 822 /* Allocate Memory for Hash Table */ 823 ierr = PetscCalloc2(ht_size,&baij->hd,ht_size,&baij->ht);CHKERRQ(ierr); 824 HD = baij->hd; 825 HT = baij->ht; 826 827 /* Loop Over A */ 828 for (i=0; i<a->mbs; i++) { 829 for (j=ai[i]; j<ai[i+1]; j++) { 830 row = i+rstart; 831 col = aj[j]+cstart; 832 833 key = row*Nbs + col + 1; 834 h1 = HASH(ht_size,key,tmp); 835 for (k=0; k<ht_size; k++) { 836 if (!HT[(h1+k)%ht_size]) { 837 HT[(h1+k)%ht_size] = key; 838 HD[(h1+k)%ht_size] = a->a + j*bs2; 839 break; 840 #if defined(PETSC_USE_INFO) 841 } else { 842 ct++; 843 #endif 844 } 845 } 846 #if defined(PETSC_USE_INFO) 847 if (k> max) max = k; 848 #endif 849 } 850 } 851 /* Loop Over B */ 852 for (i=0; i<b->mbs; i++) { 853 for (j=bi[i]; j<bi[i+1]; j++) { 854 row = i+rstart; 855 col = garray[bj[j]]; 856 key = row*Nbs + col + 1; 857 h1 = HASH(ht_size,key,tmp); 858 for (k=0; k<ht_size; k++) { 859 if (!HT[(h1+k)%ht_size]) { 860 HT[(h1+k)%ht_size] = key; 861 HD[(h1+k)%ht_size] = b->a + j*bs2; 862 break; 863 #if defined(PETSC_USE_INFO) 864 } else { 865 ct++; 866 #endif 867 } 868 } 869 #if defined(PETSC_USE_INFO) 870 if (k> max) max = k; 871 #endif 872 } 873 } 874 875 /* Print Summary */ 876 #if defined(PETSC_USE_INFO) 877 for (i=0,j=0; i<ht_size; i++) { 878 if (HT[i]) j++; 879 } 880 ierr = PetscInfo2(mat,"Average Search = %5.2f,max search = %D\n",(!j)? 0.0:((PetscReal)(ct+j))/j,max);CHKERRQ(ierr); 881 #endif 882 PetscFunctionReturn(0); 883 } 884 885 #undef __FUNCT__ 886 #define __FUNCT__ "MatAssemblyBegin_MPIBAIJ" 887 PetscErrorCode MatAssemblyBegin_MPIBAIJ(Mat mat,MatAssemblyType mode) 888 { 889 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 890 PetscErrorCode ierr; 891 PetscInt nstash,reallocs; 892 893 PetscFunctionBegin; 894 if (baij->donotstash || mat->nooffprocentries) PetscFunctionReturn(0); 895 896 ierr = MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);CHKERRQ(ierr); 897 ierr = MatStashScatterBegin_Private(mat,&mat->bstash,baij->rangebs);CHKERRQ(ierr); 898 ierr = MatStashGetInfo_Private(&mat->stash,&nstash,&reallocs);CHKERRQ(ierr); 899 ierr = PetscInfo2(mat,"Stash has %D entries,uses %D mallocs.\n",nstash,reallocs);CHKERRQ(ierr); 900 ierr = MatStashGetInfo_Private(&mat->bstash,&nstash,&reallocs);CHKERRQ(ierr); 901 ierr = PetscInfo2(mat,"Block-Stash has %D entries, uses %D mallocs.\n",nstash,reallocs);CHKERRQ(ierr); 902 PetscFunctionReturn(0); 903 } 904 905 #undef __FUNCT__ 906 #define __FUNCT__ "MatAssemblyEnd_MPIBAIJ" 907 PetscErrorCode MatAssemblyEnd_MPIBAIJ(Mat mat,MatAssemblyType mode) 908 { 909 Mat_MPIBAIJ *baij=(Mat_MPIBAIJ*)mat->data; 910 Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)baij->A->data; 911 PetscErrorCode ierr; 912 PetscInt i,j,rstart,ncols,flg,bs2=baij->bs2; 913 PetscInt *row,*col; 914 PetscBool r1,r2,r3,other_disassembled; 915 MatScalar *val; 916 PetscMPIInt n; 917 918 PetscFunctionBegin; 919 /* do not use 'b=(Mat_SeqBAIJ*)baij->B->data' as B can be reset in disassembly */ 920 if (!baij->donotstash && !mat->nooffprocentries) { 921 while (1) { 922 ierr = MatStashScatterGetMesg_Private(&mat->stash,&n,&row,&col,&val,&flg);CHKERRQ(ierr); 923 if (!flg) break; 924 925 for (i=0; i<n;) { 926 /* Now identify the consecutive vals belonging to the same row */ 927 for (j=i,rstart=row[j]; j<n; j++) { 928 if (row[j] != rstart) break; 929 } 930 if (j < n) ncols = j-i; 931 else ncols = n-i; 932 /* Now assemble all these values with a single function call */ 933 ierr = MatSetValues_MPIBAIJ(mat,1,row+i,ncols,col+i,val+i,mat->insertmode);CHKERRQ(ierr); 934 i = j; 935 } 936 } 937 ierr = MatStashScatterEnd_Private(&mat->stash);CHKERRQ(ierr); 938 /* Now process the block-stash. Since the values are stashed column-oriented, 939 set the roworiented flag to column oriented, and after MatSetValues() 940 restore the original flags */ 941 r1 = baij->roworiented; 942 r2 = a->roworiented; 943 r3 = ((Mat_SeqBAIJ*)baij->B->data)->roworiented; 944 945 baij->roworiented = PETSC_FALSE; 946 a->roworiented = PETSC_FALSE; 947 948 (((Mat_SeqBAIJ*)baij->B->data))->roworiented = PETSC_FALSE; /* b->roworiented */ 949 while (1) { 950 ierr = MatStashScatterGetMesg_Private(&mat->bstash,&n,&row,&col,&val,&flg);CHKERRQ(ierr); 951 if (!flg) break; 952 953 for (i=0; i<n;) { 954 /* Now identify the consecutive vals belonging to the same row */ 955 for (j=i,rstart=row[j]; j<n; j++) { 956 if (row[j] != rstart) break; 957 } 958 if (j < n) ncols = j-i; 959 else ncols = n-i; 960 ierr = MatSetValuesBlocked_MPIBAIJ(mat,1,row+i,ncols,col+i,val+i*bs2,mat->insertmode);CHKERRQ(ierr); 961 i = j; 962 } 963 } 964 ierr = MatStashScatterEnd_Private(&mat->bstash);CHKERRQ(ierr); 965 966 baij->roworiented = r1; 967 a->roworiented = r2; 968 969 ((Mat_SeqBAIJ*)baij->B->data)->roworiented = r3; /* b->roworiented */ 970 } 971 972 ierr = MatAssemblyBegin(baij->A,mode);CHKERRQ(ierr); 973 ierr = MatAssemblyEnd(baij->A,mode);CHKERRQ(ierr); 974 975 /* determine if any processor has disassembled, if so we must 976 also disassemble ourselfs, in order that we may reassemble. */ 977 /* 978 if nonzero structure of submatrix B cannot change then we know that 979 no processor disassembled thus we can skip this stuff 980 */ 981 if (!((Mat_SeqBAIJ*)baij->B->data)->nonew) { 982 ierr = MPIU_Allreduce(&mat->was_assembled,&other_disassembled,1,MPIU_BOOL,MPI_PROD,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 983 if (mat->was_assembled && !other_disassembled) { 984 ierr = MatDisAssemble_MPIBAIJ(mat);CHKERRQ(ierr); 985 } 986 } 987 988 if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) { 989 ierr = MatSetUpMultiply_MPIBAIJ(mat);CHKERRQ(ierr); 990 } 991 ierr = MatAssemblyBegin(baij->B,mode);CHKERRQ(ierr); 992 ierr = MatAssemblyEnd(baij->B,mode);CHKERRQ(ierr); 993 994 #if defined(PETSC_USE_INFO) 995 if (baij->ht && mode== MAT_FINAL_ASSEMBLY) { 996 ierr = PetscInfo1(mat,"Average Hash Table Search in MatSetValues = %5.2f\n",(double)((PetscReal)baij->ht_total_ct)/baij->ht_insert_ct);CHKERRQ(ierr); 997 998 baij->ht_total_ct = 0; 999 baij->ht_insert_ct = 0; 1000 } 1001 #endif 1002 if (baij->ht_flag && !baij->ht && mode == MAT_FINAL_ASSEMBLY) { 1003 ierr = MatCreateHashTable_MPIBAIJ_Private(mat,baij->ht_fact);CHKERRQ(ierr); 1004 1005 mat->ops->setvalues = MatSetValues_MPIBAIJ_HT; 1006 mat->ops->setvaluesblocked = MatSetValuesBlocked_MPIBAIJ_HT; 1007 } 1008 1009 ierr = PetscFree2(baij->rowvalues,baij->rowindices);CHKERRQ(ierr); 1010 1011 baij->rowvalues = 0; 1012 1013 /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */ 1014 if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqBAIJ*)(baij->A->data))->nonew) { 1015 PetscObjectState state = baij->A->nonzerostate + baij->B->nonzerostate; 1016 ierr = MPIU_Allreduce(&state,&mat->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1017 } 1018 PetscFunctionReturn(0); 1019 } 1020 1021 extern PetscErrorCode MatView_SeqBAIJ(Mat,PetscViewer); 1022 #include <petscdraw.h> 1023 #undef __FUNCT__ 1024 #define __FUNCT__ "MatView_MPIBAIJ_ASCIIorDraworSocket" 1025 static PetscErrorCode MatView_MPIBAIJ_ASCIIorDraworSocket(Mat mat,PetscViewer viewer) 1026 { 1027 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 1028 PetscErrorCode ierr; 1029 PetscMPIInt rank = baij->rank; 1030 PetscInt bs = mat->rmap->bs; 1031 PetscBool iascii,isdraw; 1032 PetscViewer sviewer; 1033 PetscViewerFormat format; 1034 1035 PetscFunctionBegin; 1036 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1037 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1038 if (iascii) { 1039 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1040 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1041 MatInfo info; 1042 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)mat),&rank);CHKERRQ(ierr); 1043 ierr = MatGetInfo(mat,MAT_LOCAL,&info);CHKERRQ(ierr); 1044 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr); 1045 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %D nz %D nz alloced %D bs %D mem %D\n", 1046 rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,mat->rmap->bs,(PetscInt)info.memory);CHKERRQ(ierr); 1047 ierr = MatGetInfo(baij->A,MAT_LOCAL,&info);CHKERRQ(ierr); 1048 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"[%d] on-diagonal part: nz %D \n",rank,(PetscInt)info.nz_used);CHKERRQ(ierr); 1049 ierr = MatGetInfo(baij->B,MAT_LOCAL,&info);CHKERRQ(ierr); 1050 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"[%d] off-diagonal part: nz %D \n",rank,(PetscInt)info.nz_used);CHKERRQ(ierr); 1051 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1052 ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr); 1053 ierr = PetscViewerASCIIPrintf(viewer,"Information on VecScatter used in matrix-vector product: \n");CHKERRQ(ierr); 1054 ierr = VecScatterView(baij->Mvctx,viewer);CHKERRQ(ierr); 1055 PetscFunctionReturn(0); 1056 } else if (format == PETSC_VIEWER_ASCII_INFO) { 1057 ierr = PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);CHKERRQ(ierr); 1058 PetscFunctionReturn(0); 1059 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 1060 PetscFunctionReturn(0); 1061 } 1062 } 1063 1064 if (isdraw) { 1065 PetscDraw draw; 1066 PetscBool isnull; 1067 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1068 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1069 if (isnull) PetscFunctionReturn(0); 1070 } 1071 1072 { 1073 /* assemble the entire matrix onto first processor. */ 1074 Mat A; 1075 Mat_SeqBAIJ *Aloc; 1076 PetscInt M = mat->rmap->N,N = mat->cmap->N,*ai,*aj,col,i,j,k,*rvals,mbs = baij->mbs; 1077 MatScalar *a; 1078 const char *matname; 1079 1080 /* Here we are creating a temporary matrix, so will assume MPIBAIJ is acceptable */ 1081 /* Perhaps this should be the type of mat? */ 1082 ierr = MatCreate(PetscObjectComm((PetscObject)mat),&A);CHKERRQ(ierr); 1083 if (!rank) { 1084 ierr = MatSetSizes(A,M,N,M,N);CHKERRQ(ierr); 1085 } else { 1086 ierr = MatSetSizes(A,0,0,M,N);CHKERRQ(ierr); 1087 } 1088 ierr = MatSetType(A,MATMPIBAIJ);CHKERRQ(ierr); 1089 ierr = MatMPIBAIJSetPreallocation(A,mat->rmap->bs,0,NULL,0,NULL);CHKERRQ(ierr); 1090 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 1091 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)A);CHKERRQ(ierr); 1092 1093 /* copy over the A part */ 1094 Aloc = (Mat_SeqBAIJ*)baij->A->data; 1095 ai = Aloc->i; aj = Aloc->j; a = Aloc->a; 1096 ierr = PetscMalloc1(bs,&rvals);CHKERRQ(ierr); 1097 1098 for (i=0; i<mbs; i++) { 1099 rvals[0] = bs*(baij->rstartbs + i); 1100 for (j=1; j<bs; j++) rvals[j] = rvals[j-1] + 1; 1101 for (j=ai[i]; j<ai[i+1]; j++) { 1102 col = (baij->cstartbs+aj[j])*bs; 1103 for (k=0; k<bs; k++) { 1104 ierr = MatSetValues_MPIBAIJ(A,bs,rvals,1,&col,a,INSERT_VALUES);CHKERRQ(ierr); 1105 col++; a += bs; 1106 } 1107 } 1108 } 1109 /* copy over the B part */ 1110 Aloc = (Mat_SeqBAIJ*)baij->B->data; 1111 ai = Aloc->i; aj = Aloc->j; a = Aloc->a; 1112 for (i=0; i<mbs; i++) { 1113 rvals[0] = bs*(baij->rstartbs + i); 1114 for (j=1; j<bs; j++) rvals[j] = rvals[j-1] + 1; 1115 for (j=ai[i]; j<ai[i+1]; j++) { 1116 col = baij->garray[aj[j]]*bs; 1117 for (k=0; k<bs; k++) { 1118 ierr = MatSetValues_MPIBAIJ(A,bs,rvals,1,&col,a,INSERT_VALUES);CHKERRQ(ierr); 1119 col++; a += bs; 1120 } 1121 } 1122 } 1123 ierr = PetscFree(rvals);CHKERRQ(ierr); 1124 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1125 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1126 /* 1127 Everyone has to call to draw the matrix since the graphics waits are 1128 synchronized across all processors that share the PetscDraw object 1129 */ 1130 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1131 ierr = PetscObjectGetName((PetscObject)mat,&matname);CHKERRQ(ierr); 1132 if (!rank) { 1133 ierr = PetscObjectSetName((PetscObject)((Mat_MPIBAIJ*)(A->data))->A,matname);CHKERRQ(ierr); 1134 ierr = MatView_SeqBAIJ(((Mat_MPIBAIJ*)(A->data))->A,sviewer);CHKERRQ(ierr); 1135 } 1136 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1137 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1138 ierr = MatDestroy(&A);CHKERRQ(ierr); 1139 } 1140 PetscFunctionReturn(0); 1141 } 1142 1143 #undef __FUNCT__ 1144 #define __FUNCT__ "MatView_MPIBAIJ_Binary" 1145 static PetscErrorCode MatView_MPIBAIJ_Binary(Mat mat,PetscViewer viewer) 1146 { 1147 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)mat->data; 1148 Mat_SeqBAIJ *A = (Mat_SeqBAIJ*)a->A->data; 1149 Mat_SeqBAIJ *B = (Mat_SeqBAIJ*)a->B->data; 1150 PetscErrorCode ierr; 1151 PetscInt i,*row_lens,*crow_lens,bs = mat->rmap->bs,j,k,bs2=a->bs2,header[4],nz,rlen; 1152 PetscInt *range=0,nzmax,*column_indices,cnt,col,*garray = a->garray,cstart = mat->cmap->rstart/bs,len,pcnt,l,ll; 1153 int fd; 1154 PetscScalar *column_values; 1155 FILE *file; 1156 PetscMPIInt rank,size,tag = ((PetscObject)viewer)->tag; 1157 PetscInt message_count,flowcontrolcount; 1158 1159 PetscFunctionBegin; 1160 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)mat),&rank);CHKERRQ(ierr); 1161 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRQ(ierr); 1162 nz = bs2*(A->nz + B->nz); 1163 rlen = mat->rmap->n; 1164 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1165 if (!rank) { 1166 header[0] = MAT_FILE_CLASSID; 1167 header[1] = mat->rmap->N; 1168 header[2] = mat->cmap->N; 1169 1170 ierr = MPI_Reduce(&nz,&header[3],1,MPIU_INT,MPI_SUM,0,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1171 ierr = PetscBinaryWrite(fd,header,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1172 /* get largest number of rows any processor has */ 1173 range = mat->rmap->range; 1174 for (i=1; i<size; i++) { 1175 rlen = PetscMax(rlen,range[i+1] - range[i]); 1176 } 1177 } else { 1178 ierr = MPI_Reduce(&nz,0,1,MPIU_INT,MPI_SUM,0,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1179 } 1180 1181 ierr = PetscMalloc1(rlen/bs,&crow_lens);CHKERRQ(ierr); 1182 /* compute lengths of each row */ 1183 for (i=0; i<a->mbs; i++) { 1184 crow_lens[i] = A->i[i+1] - A->i[i] + B->i[i+1] - B->i[i]; 1185 } 1186 /* store the row lengths to the file */ 1187 ierr = PetscViewerFlowControlStart(viewer,&message_count,&flowcontrolcount);CHKERRQ(ierr); 1188 if (!rank) { 1189 MPI_Status status; 1190 ierr = PetscMalloc1(rlen,&row_lens);CHKERRQ(ierr); 1191 rlen = (range[1] - range[0])/bs; 1192 for (i=0; i<rlen; i++) { 1193 for (j=0; j<bs; j++) { 1194 row_lens[i*bs+j] = bs*crow_lens[i]; 1195 } 1196 } 1197 ierr = PetscBinaryWrite(fd,row_lens,bs*rlen,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1198 for (i=1; i<size; i++) { 1199 rlen = (range[i+1] - range[i])/bs; 1200 ierr = PetscViewerFlowControlStepMaster(viewer,i,&message_count,flowcontrolcount);CHKERRQ(ierr); 1201 ierr = MPI_Recv(crow_lens,rlen,MPIU_INT,i,tag,PetscObjectComm((PetscObject)mat),&status);CHKERRQ(ierr); 1202 for (k=0; k<rlen; k++) { 1203 for (j=0; j<bs; j++) { 1204 row_lens[k*bs+j] = bs*crow_lens[k]; 1205 } 1206 } 1207 ierr = PetscBinaryWrite(fd,row_lens,bs*rlen,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1208 } 1209 ierr = PetscViewerFlowControlEndMaster(viewer,&message_count);CHKERRQ(ierr); 1210 ierr = PetscFree(row_lens);CHKERRQ(ierr); 1211 } else { 1212 ierr = PetscViewerFlowControlStepWorker(viewer,rank,&message_count);CHKERRQ(ierr); 1213 ierr = MPI_Send(crow_lens,mat->rmap->n/bs,MPIU_INT,0,tag,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1214 ierr = PetscViewerFlowControlEndWorker(viewer,&message_count);CHKERRQ(ierr); 1215 } 1216 ierr = PetscFree(crow_lens);CHKERRQ(ierr); 1217 1218 /* load up the local column indices. Include for all rows not just one for each block row since process 0 does not have the 1219 information needed to make it for each row from a block row. This does require more communication but still not more than 1220 the communication needed for the nonzero values */ 1221 nzmax = nz; /* space a largest processor needs */ 1222 ierr = MPI_Reduce(&nz,&nzmax,1,MPIU_INT,MPI_MAX,0,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1223 ierr = PetscMalloc1(nzmax,&column_indices);CHKERRQ(ierr); 1224 cnt = 0; 1225 for (i=0; i<a->mbs; i++) { 1226 pcnt = cnt; 1227 for (j=B->i[i]; j<B->i[i+1]; j++) { 1228 if ((col = garray[B->j[j]]) > cstart) break; 1229 for (l=0; l<bs; l++) { 1230 column_indices[cnt++] = bs*col+l; 1231 } 1232 } 1233 for (k=A->i[i]; k<A->i[i+1]; k++) { 1234 for (l=0; l<bs; l++) { 1235 column_indices[cnt++] = bs*(A->j[k] + cstart)+l; 1236 } 1237 } 1238 for (; j<B->i[i+1]; j++) { 1239 for (l=0; l<bs; l++) { 1240 column_indices[cnt++] = bs*garray[B->j[j]]+l; 1241 } 1242 } 1243 len = cnt - pcnt; 1244 for (k=1; k<bs; k++) { 1245 ierr = PetscMemcpy(&column_indices[cnt],&column_indices[pcnt],len*sizeof(PetscInt));CHKERRQ(ierr); 1246 cnt += len; 1247 } 1248 } 1249 if (cnt != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Internal PETSc error: cnt = %D nz = %D",cnt,nz); 1250 1251 /* store the columns to the file */ 1252 ierr = PetscViewerFlowControlStart(viewer,&message_count,&flowcontrolcount);CHKERRQ(ierr); 1253 if (!rank) { 1254 MPI_Status status; 1255 ierr = PetscBinaryWrite(fd,column_indices,nz,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1256 for (i=1; i<size; i++) { 1257 ierr = PetscViewerFlowControlStepMaster(viewer,i,&message_count,flowcontrolcount);CHKERRQ(ierr); 1258 ierr = MPI_Recv(&cnt,1,MPIU_INT,i,tag,PetscObjectComm((PetscObject)mat),&status);CHKERRQ(ierr); 1259 ierr = MPI_Recv(column_indices,cnt,MPIU_INT,i,tag,PetscObjectComm((PetscObject)mat),&status);CHKERRQ(ierr); 1260 ierr = PetscBinaryWrite(fd,column_indices,cnt,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1261 } 1262 ierr = PetscViewerFlowControlEndMaster(viewer,&message_count);CHKERRQ(ierr); 1263 } else { 1264 ierr = PetscViewerFlowControlStepWorker(viewer,rank,&message_count);CHKERRQ(ierr); 1265 ierr = MPI_Send(&cnt,1,MPIU_INT,0,tag,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1266 ierr = MPI_Send(column_indices,cnt,MPIU_INT,0,tag,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1267 ierr = PetscViewerFlowControlEndWorker(viewer,&message_count);CHKERRQ(ierr); 1268 } 1269 ierr = PetscFree(column_indices);CHKERRQ(ierr); 1270 1271 /* load up the numerical values */ 1272 ierr = PetscMalloc1(nzmax,&column_values);CHKERRQ(ierr); 1273 cnt = 0; 1274 for (i=0; i<a->mbs; i++) { 1275 rlen = bs*(B->i[i+1] - B->i[i] + A->i[i+1] - A->i[i]); 1276 for (j=B->i[i]; j<B->i[i+1]; j++) { 1277 if (garray[B->j[j]] > cstart) break; 1278 for (l=0; l<bs; l++) { 1279 for (ll=0; ll<bs; ll++) { 1280 column_values[cnt + l*rlen + ll] = B->a[bs2*j+l+bs*ll]; 1281 } 1282 } 1283 cnt += bs; 1284 } 1285 for (k=A->i[i]; k<A->i[i+1]; k++) { 1286 for (l=0; l<bs; l++) { 1287 for (ll=0; ll<bs; ll++) { 1288 column_values[cnt + l*rlen + ll] = A->a[bs2*k+l+bs*ll]; 1289 } 1290 } 1291 cnt += bs; 1292 } 1293 for (; j<B->i[i+1]; j++) { 1294 for (l=0; l<bs; l++) { 1295 for (ll=0; ll<bs; ll++) { 1296 column_values[cnt + l*rlen + ll] = B->a[bs2*j+l+bs*ll]; 1297 } 1298 } 1299 cnt += bs; 1300 } 1301 cnt += (bs-1)*rlen; 1302 } 1303 if (cnt != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal PETSc error: cnt = %D nz = %D",cnt,nz); 1304 1305 /* store the column values to the file */ 1306 ierr = PetscViewerFlowControlStart(viewer,&message_count,&flowcontrolcount);CHKERRQ(ierr); 1307 if (!rank) { 1308 MPI_Status status; 1309 ierr = PetscBinaryWrite(fd,column_values,nz,PETSC_SCALAR,PETSC_TRUE);CHKERRQ(ierr); 1310 for (i=1; i<size; i++) { 1311 ierr = PetscViewerFlowControlStepMaster(viewer,i,&message_count,flowcontrolcount);CHKERRQ(ierr); 1312 ierr = MPI_Recv(&cnt,1,MPIU_INT,i,tag,PetscObjectComm((PetscObject)mat),&status);CHKERRQ(ierr); 1313 ierr = MPI_Recv(column_values,cnt,MPIU_SCALAR,i,tag,PetscObjectComm((PetscObject)mat),&status);CHKERRQ(ierr); 1314 ierr = PetscBinaryWrite(fd,column_values,cnt,PETSC_SCALAR,PETSC_TRUE);CHKERRQ(ierr); 1315 } 1316 ierr = PetscViewerFlowControlEndMaster(viewer,&message_count);CHKERRQ(ierr); 1317 } else { 1318 ierr = PetscViewerFlowControlStepWorker(viewer,rank,&message_count);CHKERRQ(ierr); 1319 ierr = MPI_Send(&nz,1,MPIU_INT,0,tag,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1320 ierr = MPI_Send(column_values,nz,MPIU_SCALAR,0,tag,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1321 ierr = PetscViewerFlowControlEndWorker(viewer,&message_count);CHKERRQ(ierr); 1322 } 1323 ierr = PetscFree(column_values);CHKERRQ(ierr); 1324 1325 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 1326 if (file) { 1327 fprintf(file,"-matload_block_size %d\n",(int)mat->rmap->bs); 1328 } 1329 PetscFunctionReturn(0); 1330 } 1331 1332 #undef __FUNCT__ 1333 #define __FUNCT__ "MatView_MPIBAIJ" 1334 PetscErrorCode MatView_MPIBAIJ(Mat mat,PetscViewer viewer) 1335 { 1336 PetscErrorCode ierr; 1337 PetscBool iascii,isdraw,issocket,isbinary; 1338 1339 PetscFunctionBegin; 1340 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1341 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1342 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);CHKERRQ(ierr); 1343 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1344 if (iascii || isdraw || issocket) { 1345 ierr = MatView_MPIBAIJ_ASCIIorDraworSocket(mat,viewer);CHKERRQ(ierr); 1346 } else if (isbinary) { 1347 ierr = MatView_MPIBAIJ_Binary(mat,viewer);CHKERRQ(ierr); 1348 } 1349 PetscFunctionReturn(0); 1350 } 1351 1352 #undef __FUNCT__ 1353 #define __FUNCT__ "MatDestroy_MPIBAIJ" 1354 PetscErrorCode MatDestroy_MPIBAIJ(Mat mat) 1355 { 1356 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 1357 PetscErrorCode ierr; 1358 1359 PetscFunctionBegin; 1360 #if defined(PETSC_USE_LOG) 1361 PetscLogObjectState((PetscObject)mat,"Rows=%D,Cols=%D",mat->rmap->N,mat->cmap->N); 1362 #endif 1363 ierr = MatStashDestroy_Private(&mat->stash);CHKERRQ(ierr); 1364 ierr = MatStashDestroy_Private(&mat->bstash);CHKERRQ(ierr); 1365 ierr = MatDestroy(&baij->A);CHKERRQ(ierr); 1366 ierr = MatDestroy(&baij->B);CHKERRQ(ierr); 1367 #if defined(PETSC_USE_CTABLE) 1368 ierr = PetscTableDestroy(&baij->colmap);CHKERRQ(ierr); 1369 #else 1370 ierr = PetscFree(baij->colmap);CHKERRQ(ierr); 1371 #endif 1372 ierr = PetscFree(baij->garray);CHKERRQ(ierr); 1373 ierr = VecDestroy(&baij->lvec);CHKERRQ(ierr); 1374 ierr = VecScatterDestroy(&baij->Mvctx);CHKERRQ(ierr); 1375 ierr = PetscFree2(baij->rowvalues,baij->rowindices);CHKERRQ(ierr); 1376 ierr = PetscFree(baij->barray);CHKERRQ(ierr); 1377 ierr = PetscFree2(baij->hd,baij->ht);CHKERRQ(ierr); 1378 ierr = PetscFree(baij->rangebs);CHKERRQ(ierr); 1379 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1380 1381 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1382 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatStoreValues_C",NULL);CHKERRQ(ierr); 1383 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatRetrieveValues_C",NULL);CHKERRQ(ierr); 1384 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMPIBAIJSetPreallocation_C",NULL);CHKERRQ(ierr); 1385 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMPIBAIJSetPreallocationCSR_C",NULL);CHKERRQ(ierr); 1386 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",NULL);CHKERRQ(ierr); 1387 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSetHashTableFactor_C",NULL);CHKERRQ(ierr); 1388 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_mpibaij_mpisbaij_C",NULL);CHKERRQ(ierr); 1389 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_mpibaij_mpibstrm_C",NULL);CHKERRQ(ierr); 1390 PetscFunctionReturn(0); 1391 } 1392 1393 #undef __FUNCT__ 1394 #define __FUNCT__ "MatMult_MPIBAIJ" 1395 PetscErrorCode MatMult_MPIBAIJ(Mat A,Vec xx,Vec yy) 1396 { 1397 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 1398 PetscErrorCode ierr; 1399 PetscInt nt; 1400 1401 PetscFunctionBegin; 1402 ierr = VecGetLocalSize(xx,&nt);CHKERRQ(ierr); 1403 if (nt != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible partition of A and xx"); 1404 ierr = VecGetLocalSize(yy,&nt);CHKERRQ(ierr); 1405 if (nt != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible parition of A and yy"); 1406 ierr = VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1407 ierr = (*a->A->ops->mult)(a->A,xx,yy);CHKERRQ(ierr); 1408 ierr = VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1409 ierr = (*a->B->ops->multadd)(a->B,a->lvec,yy,yy);CHKERRQ(ierr); 1410 PetscFunctionReturn(0); 1411 } 1412 1413 #undef __FUNCT__ 1414 #define __FUNCT__ "MatMultAdd_MPIBAIJ" 1415 PetscErrorCode MatMultAdd_MPIBAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1416 { 1417 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 1418 PetscErrorCode ierr; 1419 1420 PetscFunctionBegin; 1421 ierr = VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1422 ierr = (*a->A->ops->multadd)(a->A,xx,yy,zz);CHKERRQ(ierr); 1423 ierr = VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1424 ierr = (*a->B->ops->multadd)(a->B,a->lvec,zz,zz);CHKERRQ(ierr); 1425 PetscFunctionReturn(0); 1426 } 1427 1428 #undef __FUNCT__ 1429 #define __FUNCT__ "MatMultTranspose_MPIBAIJ" 1430 PetscErrorCode MatMultTranspose_MPIBAIJ(Mat A,Vec xx,Vec yy) 1431 { 1432 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 1433 PetscErrorCode ierr; 1434 PetscBool merged; 1435 1436 PetscFunctionBegin; 1437 ierr = VecScatterGetMerged(a->Mvctx,&merged);CHKERRQ(ierr); 1438 /* do nondiagonal part */ 1439 ierr = (*a->B->ops->multtranspose)(a->B,xx,a->lvec);CHKERRQ(ierr); 1440 if (!merged) { 1441 /* send it on its way */ 1442 ierr = VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1443 /* do local part */ 1444 ierr = (*a->A->ops->multtranspose)(a->A,xx,yy);CHKERRQ(ierr); 1445 /* receive remote parts: note this assumes the values are not actually */ 1446 /* inserted in yy until the next line */ 1447 ierr = VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1448 } else { 1449 /* do local part */ 1450 ierr = (*a->A->ops->multtranspose)(a->A,xx,yy);CHKERRQ(ierr); 1451 /* send it on its way */ 1452 ierr = VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1453 /* values actually were received in the Begin() but we need to call this nop */ 1454 ierr = VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1455 } 1456 PetscFunctionReturn(0); 1457 } 1458 1459 #undef __FUNCT__ 1460 #define __FUNCT__ "MatMultTransposeAdd_MPIBAIJ" 1461 PetscErrorCode MatMultTransposeAdd_MPIBAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1462 { 1463 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 1464 PetscErrorCode ierr; 1465 1466 PetscFunctionBegin; 1467 /* do nondiagonal part */ 1468 ierr = (*a->B->ops->multtranspose)(a->B,xx,a->lvec);CHKERRQ(ierr); 1469 /* send it on its way */ 1470 ierr = VecScatterBegin(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1471 /* do local part */ 1472 ierr = (*a->A->ops->multtransposeadd)(a->A,xx,yy,zz);CHKERRQ(ierr); 1473 /* receive remote parts: note this assumes the values are not actually */ 1474 /* inserted in yy until the next line, which is true for my implementation*/ 1475 /* but is not perhaps always true. */ 1476 ierr = VecScatterEnd(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1477 PetscFunctionReturn(0); 1478 } 1479 1480 /* 1481 This only works correctly for square matrices where the subblock A->A is the 1482 diagonal block 1483 */ 1484 #undef __FUNCT__ 1485 #define __FUNCT__ "MatGetDiagonal_MPIBAIJ" 1486 PetscErrorCode MatGetDiagonal_MPIBAIJ(Mat A,Vec v) 1487 { 1488 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 1489 PetscErrorCode ierr; 1490 1491 PetscFunctionBegin; 1492 if (A->rmap->N != A->cmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Supports only square matrix where A->A is diag block"); 1493 ierr = MatGetDiagonal(a->A,v);CHKERRQ(ierr); 1494 PetscFunctionReturn(0); 1495 } 1496 1497 #undef __FUNCT__ 1498 #define __FUNCT__ "MatScale_MPIBAIJ" 1499 PetscErrorCode MatScale_MPIBAIJ(Mat A,PetscScalar aa) 1500 { 1501 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 1502 PetscErrorCode ierr; 1503 1504 PetscFunctionBegin; 1505 ierr = MatScale(a->A,aa);CHKERRQ(ierr); 1506 ierr = MatScale(a->B,aa);CHKERRQ(ierr); 1507 PetscFunctionReturn(0); 1508 } 1509 1510 #undef __FUNCT__ 1511 #define __FUNCT__ "MatGetRow_MPIBAIJ" 1512 PetscErrorCode MatGetRow_MPIBAIJ(Mat matin,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1513 { 1514 Mat_MPIBAIJ *mat = (Mat_MPIBAIJ*)matin->data; 1515 PetscScalar *vworkA,*vworkB,**pvA,**pvB,*v_p; 1516 PetscErrorCode ierr; 1517 PetscInt bs = matin->rmap->bs,bs2 = mat->bs2,i,*cworkA,*cworkB,**pcA,**pcB; 1518 PetscInt nztot,nzA,nzB,lrow,brstart = matin->rmap->rstart,brend = matin->rmap->rend; 1519 PetscInt *cmap,*idx_p,cstart = mat->cstartbs; 1520 1521 PetscFunctionBegin; 1522 if (row < brstart || row >= brend) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only local rows"); 1523 if (mat->getrowactive) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Already active"); 1524 mat->getrowactive = PETSC_TRUE; 1525 1526 if (!mat->rowvalues && (idx || v)) { 1527 /* 1528 allocate enough space to hold information from the longest row. 1529 */ 1530 Mat_SeqBAIJ *Aa = (Mat_SeqBAIJ*)mat->A->data,*Ba = (Mat_SeqBAIJ*)mat->B->data; 1531 PetscInt max = 1,mbs = mat->mbs,tmp; 1532 for (i=0; i<mbs; i++) { 1533 tmp = Aa->i[i+1] - Aa->i[i] + Ba->i[i+1] - Ba->i[i]; 1534 if (max < tmp) max = tmp; 1535 } 1536 ierr = PetscMalloc2(max*bs2,&mat->rowvalues,max*bs2,&mat->rowindices);CHKERRQ(ierr); 1537 } 1538 lrow = row - brstart; 1539 1540 pvA = &vworkA; pcA = &cworkA; pvB = &vworkB; pcB = &cworkB; 1541 if (!v) {pvA = 0; pvB = 0;} 1542 if (!idx) {pcA = 0; if (!v) pcB = 0;} 1543 ierr = (*mat->A->ops->getrow)(mat->A,lrow,&nzA,pcA,pvA);CHKERRQ(ierr); 1544 ierr = (*mat->B->ops->getrow)(mat->B,lrow,&nzB,pcB,pvB);CHKERRQ(ierr); 1545 nztot = nzA + nzB; 1546 1547 cmap = mat->garray; 1548 if (v || idx) { 1549 if (nztot) { 1550 /* Sort by increasing column numbers, assuming A and B already sorted */ 1551 PetscInt imark = -1; 1552 if (v) { 1553 *v = v_p = mat->rowvalues; 1554 for (i=0; i<nzB; i++) { 1555 if (cmap[cworkB[i]/bs] < cstart) v_p[i] = vworkB[i]; 1556 else break; 1557 } 1558 imark = i; 1559 for (i=0; i<nzA; i++) v_p[imark+i] = vworkA[i]; 1560 for (i=imark; i<nzB; i++) v_p[nzA+i] = vworkB[i]; 1561 } 1562 if (idx) { 1563 *idx = idx_p = mat->rowindices; 1564 if (imark > -1) { 1565 for (i=0; i<imark; i++) { 1566 idx_p[i] = cmap[cworkB[i]/bs]*bs + cworkB[i]%bs; 1567 } 1568 } else { 1569 for (i=0; i<nzB; i++) { 1570 if (cmap[cworkB[i]/bs] < cstart) idx_p[i] = cmap[cworkB[i]/bs]*bs + cworkB[i]%bs; 1571 else break; 1572 } 1573 imark = i; 1574 } 1575 for (i=0; i<nzA; i++) idx_p[imark+i] = cstart*bs + cworkA[i]; 1576 for (i=imark; i<nzB; i++) idx_p[nzA+i] = cmap[cworkB[i]/bs]*bs + cworkB[i]%bs ; 1577 } 1578 } else { 1579 if (idx) *idx = 0; 1580 if (v) *v = 0; 1581 } 1582 } 1583 *nz = nztot; 1584 ierr = (*mat->A->ops->restorerow)(mat->A,lrow,&nzA,pcA,pvA);CHKERRQ(ierr); 1585 ierr = (*mat->B->ops->restorerow)(mat->B,lrow,&nzB,pcB,pvB);CHKERRQ(ierr); 1586 PetscFunctionReturn(0); 1587 } 1588 1589 #undef __FUNCT__ 1590 #define __FUNCT__ "MatRestoreRow_MPIBAIJ" 1591 PetscErrorCode MatRestoreRow_MPIBAIJ(Mat mat,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1592 { 1593 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 1594 1595 PetscFunctionBegin; 1596 if (!baij->getrowactive) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"MatGetRow not called"); 1597 baij->getrowactive = PETSC_FALSE; 1598 PetscFunctionReturn(0); 1599 } 1600 1601 #undef __FUNCT__ 1602 #define __FUNCT__ "MatZeroEntries_MPIBAIJ" 1603 PetscErrorCode MatZeroEntries_MPIBAIJ(Mat A) 1604 { 1605 Mat_MPIBAIJ *l = (Mat_MPIBAIJ*)A->data; 1606 PetscErrorCode ierr; 1607 1608 PetscFunctionBegin; 1609 ierr = MatZeroEntries(l->A);CHKERRQ(ierr); 1610 ierr = MatZeroEntries(l->B);CHKERRQ(ierr); 1611 PetscFunctionReturn(0); 1612 } 1613 1614 #undef __FUNCT__ 1615 #define __FUNCT__ "MatGetInfo_MPIBAIJ" 1616 PetscErrorCode MatGetInfo_MPIBAIJ(Mat matin,MatInfoType flag,MatInfo *info) 1617 { 1618 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)matin->data; 1619 Mat A = a->A,B = a->B; 1620 PetscErrorCode ierr; 1621 PetscReal isend[5],irecv[5]; 1622 1623 PetscFunctionBegin; 1624 info->block_size = (PetscReal)matin->rmap->bs; 1625 1626 ierr = MatGetInfo(A,MAT_LOCAL,info);CHKERRQ(ierr); 1627 1628 isend[0] = info->nz_used; isend[1] = info->nz_allocated; isend[2] = info->nz_unneeded; 1629 isend[3] = info->memory; isend[4] = info->mallocs; 1630 1631 ierr = MatGetInfo(B,MAT_LOCAL,info);CHKERRQ(ierr); 1632 1633 isend[0] += info->nz_used; isend[1] += info->nz_allocated; isend[2] += info->nz_unneeded; 1634 isend[3] += info->memory; isend[4] += info->mallocs; 1635 1636 if (flag == MAT_LOCAL) { 1637 info->nz_used = isend[0]; 1638 info->nz_allocated = isend[1]; 1639 info->nz_unneeded = isend[2]; 1640 info->memory = isend[3]; 1641 info->mallocs = isend[4]; 1642 } else if (flag == MAT_GLOBAL_MAX) { 1643 ierr = MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)matin));CHKERRQ(ierr); 1644 1645 info->nz_used = irecv[0]; 1646 info->nz_allocated = irecv[1]; 1647 info->nz_unneeded = irecv[2]; 1648 info->memory = irecv[3]; 1649 info->mallocs = irecv[4]; 1650 } else if (flag == MAT_GLOBAL_SUM) { 1651 ierr = MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)matin));CHKERRQ(ierr); 1652 1653 info->nz_used = irecv[0]; 1654 info->nz_allocated = irecv[1]; 1655 info->nz_unneeded = irecv[2]; 1656 info->memory = irecv[3]; 1657 info->mallocs = irecv[4]; 1658 } else SETERRQ1(PetscObjectComm((PetscObject)matin),PETSC_ERR_ARG_WRONG,"Unknown MatInfoType argument %d",(int)flag); 1659 info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */ 1660 info->fill_ratio_needed = 0; 1661 info->factor_mallocs = 0; 1662 PetscFunctionReturn(0); 1663 } 1664 1665 #undef __FUNCT__ 1666 #define __FUNCT__ "MatSetOption_MPIBAIJ" 1667 PetscErrorCode MatSetOption_MPIBAIJ(Mat A,MatOption op,PetscBool flg) 1668 { 1669 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 1670 PetscErrorCode ierr; 1671 1672 PetscFunctionBegin; 1673 switch (op) { 1674 case MAT_NEW_NONZERO_LOCATIONS: 1675 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1676 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1677 case MAT_KEEP_NONZERO_PATTERN: 1678 case MAT_NEW_NONZERO_LOCATION_ERR: 1679 MatCheckPreallocated(A,1); 1680 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 1681 ierr = MatSetOption(a->B,op,flg);CHKERRQ(ierr); 1682 break; 1683 case MAT_ROW_ORIENTED: 1684 MatCheckPreallocated(A,1); 1685 a->roworiented = flg; 1686 1687 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 1688 ierr = MatSetOption(a->B,op,flg);CHKERRQ(ierr); 1689 break; 1690 case MAT_NEW_DIAGONALS: 1691 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1692 break; 1693 case MAT_IGNORE_OFF_PROC_ENTRIES: 1694 a->donotstash = flg; 1695 break; 1696 case MAT_USE_HASH_TABLE: 1697 a->ht_flag = flg; 1698 a->ht_fact = 1.39; 1699 break; 1700 case MAT_SYMMETRIC: 1701 case MAT_STRUCTURALLY_SYMMETRIC: 1702 case MAT_HERMITIAN: 1703 case MAT_SUBMAT_SINGLEIS: 1704 case MAT_SYMMETRY_ETERNAL: 1705 MatCheckPreallocated(A,1); 1706 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 1707 break; 1708 default: 1709 SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"unknown option %d",op); 1710 } 1711 PetscFunctionReturn(0); 1712 } 1713 1714 #undef __FUNCT__ 1715 #define __FUNCT__ "MatTranspose_MPIBAIJ" 1716 PetscErrorCode MatTranspose_MPIBAIJ(Mat A,MatReuse reuse,Mat *matout) 1717 { 1718 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)A->data; 1719 Mat_SeqBAIJ *Aloc; 1720 Mat B; 1721 PetscErrorCode ierr; 1722 PetscInt M =A->rmap->N,N=A->cmap->N,*ai,*aj,i,*rvals,j,k,col; 1723 PetscInt bs=A->rmap->bs,mbs=baij->mbs; 1724 MatScalar *a; 1725 1726 PetscFunctionBegin; 1727 if (reuse == MAT_REUSE_MATRIX && A == *matout && M != N) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1728 if (reuse == MAT_INITIAL_MATRIX || *matout == A) { 1729 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1730 ierr = MatSetSizes(B,A->cmap->n,A->rmap->n,N,M);CHKERRQ(ierr); 1731 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1732 /* Do not know preallocation information, but must set block size */ 1733 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,PETSC_DECIDE,NULL,PETSC_DECIDE,NULL);CHKERRQ(ierr); 1734 } else { 1735 B = *matout; 1736 } 1737 1738 /* copy over the A part */ 1739 Aloc = (Mat_SeqBAIJ*)baij->A->data; 1740 ai = Aloc->i; aj = Aloc->j; a = Aloc->a; 1741 ierr = PetscMalloc1(bs,&rvals);CHKERRQ(ierr); 1742 1743 for (i=0; i<mbs; i++) { 1744 rvals[0] = bs*(baij->rstartbs + i); 1745 for (j=1; j<bs; j++) rvals[j] = rvals[j-1] + 1; 1746 for (j=ai[i]; j<ai[i+1]; j++) { 1747 col = (baij->cstartbs+aj[j])*bs; 1748 for (k=0; k<bs; k++) { 1749 ierr = MatSetValues_MPIBAIJ(B,1,&col,bs,rvals,a,INSERT_VALUES);CHKERRQ(ierr); 1750 1751 col++; a += bs; 1752 } 1753 } 1754 } 1755 /* copy over the B part */ 1756 Aloc = (Mat_SeqBAIJ*)baij->B->data; 1757 ai = Aloc->i; aj = Aloc->j; a = Aloc->a; 1758 for (i=0; i<mbs; i++) { 1759 rvals[0] = bs*(baij->rstartbs + i); 1760 for (j=1; j<bs; j++) rvals[j] = rvals[j-1] + 1; 1761 for (j=ai[i]; j<ai[i+1]; j++) { 1762 col = baij->garray[aj[j]]*bs; 1763 for (k=0; k<bs; k++) { 1764 ierr = MatSetValues_MPIBAIJ(B,1,&col,bs,rvals,a,INSERT_VALUES);CHKERRQ(ierr); 1765 col++; 1766 a += bs; 1767 } 1768 } 1769 } 1770 ierr = PetscFree(rvals);CHKERRQ(ierr); 1771 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1772 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1773 1774 if (reuse == MAT_INITIAL_MATRIX || *matout != A) *matout = B; 1775 else { 1776 ierr = MatHeaderMerge(A,&B);CHKERRQ(ierr); 1777 } 1778 PetscFunctionReturn(0); 1779 } 1780 1781 #undef __FUNCT__ 1782 #define __FUNCT__ "MatDiagonalScale_MPIBAIJ" 1783 PetscErrorCode MatDiagonalScale_MPIBAIJ(Mat mat,Vec ll,Vec rr) 1784 { 1785 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 1786 Mat a = baij->A,b = baij->B; 1787 PetscErrorCode ierr; 1788 PetscInt s1,s2,s3; 1789 1790 PetscFunctionBegin; 1791 ierr = MatGetLocalSize(mat,&s2,&s3);CHKERRQ(ierr); 1792 if (rr) { 1793 ierr = VecGetLocalSize(rr,&s1);CHKERRQ(ierr); 1794 if (s1!=s3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"right vector non-conforming local size"); 1795 /* Overlap communication with computation. */ 1796 ierr = VecScatterBegin(baij->Mvctx,rr,baij->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1797 } 1798 if (ll) { 1799 ierr = VecGetLocalSize(ll,&s1);CHKERRQ(ierr); 1800 if (s1!=s2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"left vector non-conforming local size"); 1801 ierr = (*b->ops->diagonalscale)(b,ll,NULL);CHKERRQ(ierr); 1802 } 1803 /* scale the diagonal block */ 1804 ierr = (*a->ops->diagonalscale)(a,ll,rr);CHKERRQ(ierr); 1805 1806 if (rr) { 1807 /* Do a scatter end and then right scale the off-diagonal block */ 1808 ierr = VecScatterEnd(baij->Mvctx,rr,baij->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1809 ierr = (*b->ops->diagonalscale)(b,NULL,baij->lvec);CHKERRQ(ierr); 1810 } 1811 PetscFunctionReturn(0); 1812 } 1813 1814 #undef __FUNCT__ 1815 #define __FUNCT__ "MatZeroRows_MPIBAIJ" 1816 PetscErrorCode MatZeroRows_MPIBAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1817 { 1818 Mat_MPIBAIJ *l = (Mat_MPIBAIJ *) A->data; 1819 PetscInt *lrows; 1820 PetscInt r, len; 1821 PetscErrorCode ierr; 1822 1823 PetscFunctionBegin; 1824 /* get locally owned rows */ 1825 ierr = MatZeroRowsMapLocal_Private(A,N,rows,&len,&lrows);CHKERRQ(ierr); 1826 /* fix right hand side if needed */ 1827 if (x && b) { 1828 const PetscScalar *xx; 1829 PetscScalar *bb; 1830 1831 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1832 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1833 for (r = 0; r < len; ++r) bb[lrows[r]] = diag*xx[lrows[r]]; 1834 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1835 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1836 } 1837 1838 /* actually zap the local rows */ 1839 /* 1840 Zero the required rows. If the "diagonal block" of the matrix 1841 is square and the user wishes to set the diagonal we use separate 1842 code so that MatSetValues() is not called for each diagonal allocating 1843 new memory, thus calling lots of mallocs and slowing things down. 1844 1845 */ 1846 /* must zero l->B before l->A because the (diag) case below may put values into l->B*/ 1847 ierr = MatZeroRows_SeqBAIJ(l->B,len,lrows,0.0,NULL,NULL);CHKERRQ(ierr); 1848 if (A->congruentlayouts == -1) { /* first time we compare rows and cols layouts */ 1849 PetscBool cong; 1850 ierr = PetscLayoutCompare(A->rmap,A->cmap,&cong);CHKERRQ(ierr); 1851 if (cong) A->congruentlayouts = 1; 1852 else A->congruentlayouts = 0; 1853 } 1854 if ((diag != 0.0) && A->congruentlayouts) { 1855 ierr = MatZeroRows_SeqBAIJ(l->A,len,lrows,diag,NULL,NULL);CHKERRQ(ierr); 1856 } else if (diag != 0.0) { 1857 ierr = MatZeroRows_SeqBAIJ(l->A,len,lrows,0.0,0,0);CHKERRQ(ierr); 1858 if (((Mat_SeqBAIJ*)l->A->data)->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatZeroRows() on rectangular matrices cannot be used with the Mat options \n\ 1859 MAT_NEW_NONZERO_LOCATIONS,MAT_NEW_NONZERO_LOCATION_ERR,MAT_NEW_NONZERO_ALLOCATION_ERR"); 1860 for (r = 0; r < len; ++r) { 1861 const PetscInt row = lrows[r] + A->rmap->rstart; 1862 ierr = MatSetValues(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr); 1863 } 1864 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1865 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1866 } else { 1867 ierr = MatZeroRows_SeqBAIJ(l->A,len,lrows,0.0,NULL,NULL);CHKERRQ(ierr); 1868 } 1869 ierr = PetscFree(lrows);CHKERRQ(ierr); 1870 1871 /* only change matrix nonzero state if pattern was allowed to be changed */ 1872 if (!((Mat_SeqBAIJ*)(l->A->data))->keepnonzeropattern) { 1873 PetscObjectState state = l->A->nonzerostate + l->B->nonzerostate; 1874 ierr = MPIU_Allreduce(&state,&A->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1875 } 1876 PetscFunctionReturn(0); 1877 } 1878 1879 #undef __FUNCT__ 1880 #define __FUNCT__ "MatZeroRowsColumns_MPIBAIJ" 1881 PetscErrorCode MatZeroRowsColumns_MPIBAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1882 { 1883 Mat_MPIBAIJ *l = (Mat_MPIBAIJ*)A->data; 1884 PetscErrorCode ierr; 1885 PetscMPIInt n = A->rmap->n; 1886 PetscInt i,j,k,r,p = 0,len = 0,row,col,count; 1887 PetscInt *lrows,*owners = A->rmap->range; 1888 PetscSFNode *rrows; 1889 PetscSF sf; 1890 const PetscScalar *xx; 1891 PetscScalar *bb,*mask; 1892 Vec xmask,lmask; 1893 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ*)l->B->data; 1894 PetscInt bs = A->rmap->bs, bs2 = baij->bs2; 1895 PetscScalar *aa; 1896 1897 PetscFunctionBegin; 1898 /* Create SF where leaves are input rows and roots are owned rows */ 1899 ierr = PetscMalloc1(n, &lrows);CHKERRQ(ierr); 1900 for (r = 0; r < n; ++r) lrows[r] = -1; 1901 ierr = PetscMalloc1(N, &rrows);CHKERRQ(ierr); 1902 for (r = 0; r < N; ++r) { 1903 const PetscInt idx = rows[r]; 1904 if (idx < 0 || A->rmap->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range [0,%D)",idx,A->rmap->N); 1905 if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this row too */ 1906 ierr = PetscLayoutFindOwner(A->rmap,idx,&p);CHKERRQ(ierr); 1907 } 1908 rrows[r].rank = p; 1909 rrows[r].index = rows[r] - owners[p]; 1910 } 1911 ierr = PetscSFCreate(PetscObjectComm((PetscObject) A), &sf);CHKERRQ(ierr); 1912 ierr = PetscSFSetGraph(sf, n, N, NULL, PETSC_OWN_POINTER, rrows, PETSC_OWN_POINTER);CHKERRQ(ierr); 1913 /* Collect flags for rows to be zeroed */ 1914 ierr = PetscSFReduceBegin(sf, MPIU_INT, (PetscInt *) rows, lrows, MPI_LOR);CHKERRQ(ierr); 1915 ierr = PetscSFReduceEnd(sf, MPIU_INT, (PetscInt *) rows, lrows, MPI_LOR);CHKERRQ(ierr); 1916 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 1917 /* Compress and put in row numbers */ 1918 for (r = 0; r < n; ++r) if (lrows[r] >= 0) lrows[len++] = r; 1919 /* zero diagonal part of matrix */ 1920 ierr = MatZeroRowsColumns(l->A,len,lrows,diag,x,b);CHKERRQ(ierr); 1921 /* handle off diagonal part of matrix */ 1922 ierr = MatCreateVecs(A,&xmask,NULL);CHKERRQ(ierr); 1923 ierr = VecDuplicate(l->lvec,&lmask);CHKERRQ(ierr); 1924 ierr = VecGetArray(xmask,&bb);CHKERRQ(ierr); 1925 for (i=0; i<len; i++) bb[lrows[i]] = 1; 1926 ierr = VecRestoreArray(xmask,&bb);CHKERRQ(ierr); 1927 ierr = VecScatterBegin(l->Mvctx,xmask,lmask,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1928 ierr = VecScatterEnd(l->Mvctx,xmask,lmask,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1929 ierr = VecDestroy(&xmask);CHKERRQ(ierr); 1930 if (x) { 1931 ierr = VecScatterBegin(l->Mvctx,x,l->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1932 ierr = VecScatterEnd(l->Mvctx,x,l->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1933 ierr = VecGetArrayRead(l->lvec,&xx);CHKERRQ(ierr); 1934 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1935 } 1936 ierr = VecGetArray(lmask,&mask);CHKERRQ(ierr); 1937 /* remove zeroed rows of off diagonal matrix */ 1938 for (i = 0; i < len; ++i) { 1939 row = lrows[i]; 1940 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1941 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1942 for (k = 0; k < count; ++k) { 1943 aa[0] = 0.0; 1944 aa += bs; 1945 } 1946 } 1947 /* loop over all elements of off process part of matrix zeroing removed columns*/ 1948 for (i = 0; i < l->B->rmap->N; ++i) { 1949 row = i/bs; 1950 for (j = baij->i[row]; j < baij->i[row+1]; ++j) { 1951 for (k = 0; k < bs; ++k) { 1952 col = bs*baij->j[j] + k; 1953 if (PetscAbsScalar(mask[col])) { 1954 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1955 if (x) bb[i] -= aa[0]*xx[col]; 1956 aa[0] = 0.0; 1957 } 1958 } 1959 } 1960 } 1961 if (x) { 1962 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1963 ierr = VecRestoreArrayRead(l->lvec,&xx);CHKERRQ(ierr); 1964 } 1965 ierr = VecRestoreArray(lmask,&mask);CHKERRQ(ierr); 1966 ierr = VecDestroy(&lmask);CHKERRQ(ierr); 1967 ierr = PetscFree(lrows);CHKERRQ(ierr); 1968 1969 /* only change matrix nonzero state if pattern was allowed to be changed */ 1970 if (!((Mat_SeqBAIJ*)(l->A->data))->keepnonzeropattern) { 1971 PetscObjectState state = l->A->nonzerostate + l->B->nonzerostate; 1972 ierr = MPIU_Allreduce(&state,&A->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1973 } 1974 PetscFunctionReturn(0); 1975 } 1976 1977 #undef __FUNCT__ 1978 #define __FUNCT__ "MatSetUnfactored_MPIBAIJ" 1979 PetscErrorCode MatSetUnfactored_MPIBAIJ(Mat A) 1980 { 1981 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 1982 PetscErrorCode ierr; 1983 1984 PetscFunctionBegin; 1985 ierr = MatSetUnfactored(a->A);CHKERRQ(ierr); 1986 PetscFunctionReturn(0); 1987 } 1988 1989 static PetscErrorCode MatDuplicate_MPIBAIJ(Mat,MatDuplicateOption,Mat*); 1990 1991 #undef __FUNCT__ 1992 #define __FUNCT__ "MatEqual_MPIBAIJ" 1993 PetscErrorCode MatEqual_MPIBAIJ(Mat A,Mat B,PetscBool *flag) 1994 { 1995 Mat_MPIBAIJ *matB = (Mat_MPIBAIJ*)B->data,*matA = (Mat_MPIBAIJ*)A->data; 1996 Mat a,b,c,d; 1997 PetscBool flg; 1998 PetscErrorCode ierr; 1999 2000 PetscFunctionBegin; 2001 a = matA->A; b = matA->B; 2002 c = matB->A; d = matB->B; 2003 2004 ierr = MatEqual(a,c,&flg);CHKERRQ(ierr); 2005 if (flg) { 2006 ierr = MatEqual(b,d,&flg);CHKERRQ(ierr); 2007 } 2008 ierr = MPIU_Allreduce(&flg,flag,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2009 PetscFunctionReturn(0); 2010 } 2011 2012 #undef __FUNCT__ 2013 #define __FUNCT__ "MatCopy_MPIBAIJ" 2014 PetscErrorCode MatCopy_MPIBAIJ(Mat A,Mat B,MatStructure str) 2015 { 2016 PetscErrorCode ierr; 2017 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 2018 Mat_MPIBAIJ *b = (Mat_MPIBAIJ*)B->data; 2019 2020 PetscFunctionBegin; 2021 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2022 if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) { 2023 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2024 } else { 2025 ierr = MatCopy(a->A,b->A,str);CHKERRQ(ierr); 2026 ierr = MatCopy(a->B,b->B,str);CHKERRQ(ierr); 2027 } 2028 PetscFunctionReturn(0); 2029 } 2030 2031 #undef __FUNCT__ 2032 #define __FUNCT__ "MatSetUp_MPIBAIJ" 2033 PetscErrorCode MatSetUp_MPIBAIJ(Mat A) 2034 { 2035 PetscErrorCode ierr; 2036 2037 PetscFunctionBegin; 2038 ierr = MatMPIBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,0,PETSC_DEFAULT,0);CHKERRQ(ierr); 2039 PetscFunctionReturn(0); 2040 } 2041 2042 #undef __FUNCT__ 2043 #define __FUNCT__ "MatAXPYGetPreallocation_MPIBAIJ" 2044 PetscErrorCode MatAXPYGetPreallocation_MPIBAIJ(Mat Y,const PetscInt *yltog,Mat X,const PetscInt *xltog,PetscInt *nnz) 2045 { 2046 PetscErrorCode ierr; 2047 PetscInt bs = Y->rmap->bs,m = Y->rmap->N/bs; 2048 Mat_SeqBAIJ *x = (Mat_SeqBAIJ*)X->data; 2049 Mat_SeqBAIJ *y = (Mat_SeqBAIJ*)Y->data; 2050 2051 PetscFunctionBegin; 2052 ierr = MatAXPYGetPreallocation_MPIX_private(m,x->i,x->j,xltog,y->i,y->j,yltog,nnz);CHKERRQ(ierr); 2053 PetscFunctionReturn(0); 2054 } 2055 2056 #undef __FUNCT__ 2057 #define __FUNCT__ "MatAXPY_MPIBAIJ" 2058 PetscErrorCode MatAXPY_MPIBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2059 { 2060 PetscErrorCode ierr; 2061 Mat_MPIBAIJ *xx=(Mat_MPIBAIJ*)X->data,*yy=(Mat_MPIBAIJ*)Y->data; 2062 PetscBLASInt bnz,one=1; 2063 Mat_SeqBAIJ *x,*y; 2064 2065 PetscFunctionBegin; 2066 if (str == SAME_NONZERO_PATTERN) { 2067 PetscScalar alpha = a; 2068 x = (Mat_SeqBAIJ*)xx->A->data; 2069 y = (Mat_SeqBAIJ*)yy->A->data; 2070 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2071 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2072 x = (Mat_SeqBAIJ*)xx->B->data; 2073 y = (Mat_SeqBAIJ*)yy->B->data; 2074 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2075 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2076 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2077 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2078 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2079 } else { 2080 Mat B; 2081 PetscInt *nnz_d,*nnz_o,bs=Y->rmap->bs; 2082 ierr = PetscMalloc1(yy->A->rmap->N,&nnz_d);CHKERRQ(ierr); 2083 ierr = PetscMalloc1(yy->B->rmap->N,&nnz_o);CHKERRQ(ierr); 2084 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2085 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2086 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2087 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2088 ierr = MatSetType(B,MATMPIBAIJ);CHKERRQ(ierr); 2089 ierr = MatAXPYGetPreallocation_SeqBAIJ(yy->A,xx->A,nnz_d);CHKERRQ(ierr); 2090 ierr = MatAXPYGetPreallocation_MPIBAIJ(yy->B,yy->garray,xx->B,xx->garray,nnz_o);CHKERRQ(ierr); 2091 ierr = MatMPIBAIJSetPreallocation(B,bs,0,nnz_d,0,nnz_o);CHKERRQ(ierr); 2092 /* MatAXPY_BasicWithPreallocation() for BAIJ matrix is much slower than AIJ, even for bs=1 ! */ 2093 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2094 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2095 ierr = PetscFree(nnz_d);CHKERRQ(ierr); 2096 ierr = PetscFree(nnz_o);CHKERRQ(ierr); 2097 } 2098 PetscFunctionReturn(0); 2099 } 2100 2101 #undef __FUNCT__ 2102 #define __FUNCT__ "MatRealPart_MPIBAIJ" 2103 PetscErrorCode MatRealPart_MPIBAIJ(Mat A) 2104 { 2105 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 2106 PetscErrorCode ierr; 2107 2108 PetscFunctionBegin; 2109 ierr = MatRealPart(a->A);CHKERRQ(ierr); 2110 ierr = MatRealPart(a->B);CHKERRQ(ierr); 2111 PetscFunctionReturn(0); 2112 } 2113 2114 #undef __FUNCT__ 2115 #define __FUNCT__ "MatImaginaryPart_MPIBAIJ" 2116 PetscErrorCode MatImaginaryPart_MPIBAIJ(Mat A) 2117 { 2118 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 2119 PetscErrorCode ierr; 2120 2121 PetscFunctionBegin; 2122 ierr = MatImaginaryPart(a->A);CHKERRQ(ierr); 2123 ierr = MatImaginaryPart(a->B);CHKERRQ(ierr); 2124 PetscFunctionReturn(0); 2125 } 2126 2127 #undef __FUNCT__ 2128 #define __FUNCT__ "MatGetSubMatrix_MPIBAIJ" 2129 PetscErrorCode MatGetSubMatrix_MPIBAIJ(Mat mat,IS isrow,IS iscol,MatReuse call,Mat *newmat) 2130 { 2131 PetscErrorCode ierr; 2132 IS iscol_local; 2133 PetscInt csize; 2134 2135 PetscFunctionBegin; 2136 ierr = ISGetLocalSize(iscol,&csize);CHKERRQ(ierr); 2137 if (call == MAT_REUSE_MATRIX) { 2138 ierr = PetscObjectQuery((PetscObject)*newmat,"ISAllGather",(PetscObject*)&iscol_local);CHKERRQ(ierr); 2139 if (!iscol_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Submatrix passed in was not used before, cannot reuse"); 2140 } else { 2141 ierr = ISAllGather(iscol,&iscol_local);CHKERRQ(ierr); 2142 } 2143 ierr = MatGetSubMatrix_MPIBAIJ_Private(mat,isrow,iscol_local,csize,call,newmat);CHKERRQ(ierr); 2144 if (call == MAT_INITIAL_MATRIX) { 2145 ierr = PetscObjectCompose((PetscObject)*newmat,"ISAllGather",(PetscObject)iscol_local);CHKERRQ(ierr); 2146 ierr = ISDestroy(&iscol_local);CHKERRQ(ierr); 2147 } 2148 PetscFunctionReturn(0); 2149 } 2150 extern PetscErrorCode MatGetSubMatrices_MPIBAIJ_local(Mat,PetscInt,const IS[],const IS[],MatReuse,PetscBool*,PetscBool*,Mat*); 2151 #undef __FUNCT__ 2152 #define __FUNCT__ "MatGetSubMatrix_MPIBAIJ_Private" 2153 /* 2154 Not great since it makes two copies of the submatrix, first an SeqBAIJ 2155 in local and then by concatenating the local matrices the end result. 2156 Writing it directly would be much like MatGetSubMatrices_MPIBAIJ(). 2157 This routine is used for BAIJ and SBAIJ matrices (unfortunate dependency). 2158 */ 2159 PetscErrorCode MatGetSubMatrix_MPIBAIJ_Private(Mat mat,IS isrow,IS iscol,PetscInt csize,MatReuse call,Mat *newmat) 2160 { 2161 PetscErrorCode ierr; 2162 PetscMPIInt rank,size; 2163 PetscInt i,m,n,rstart,row,rend,nz,*cwork,j,bs; 2164 PetscInt *ii,*jj,nlocal,*dlens,*olens,dlen,olen,jend,mglobal,ncol,nrow; 2165 Mat M,Mreuse; 2166 MatScalar *vwork,*aa; 2167 MPI_Comm comm; 2168 IS isrow_new, iscol_new; 2169 PetscBool idflag,allrows, allcols; 2170 Mat_SeqBAIJ *aij; 2171 2172 PetscFunctionBegin; 2173 ierr = PetscObjectGetComm((PetscObject)mat,&comm);CHKERRQ(ierr); 2174 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 2175 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2176 /* The compression and expansion should be avoided. Doesn't point 2177 out errors, might change the indices, hence buggey */ 2178 ierr = ISCompressIndicesGeneral(mat->rmap->N,mat->rmap->n,mat->rmap->bs,1,&isrow,&isrow_new);CHKERRQ(ierr); 2179 ierr = ISCompressIndicesGeneral(mat->cmap->N,mat->cmap->n,mat->cmap->bs,1,&iscol,&iscol_new);CHKERRQ(ierr); 2180 2181 /* Check for special case: each processor gets entire matrix columns */ 2182 ierr = ISIdentity(iscol,&idflag);CHKERRQ(ierr); 2183 ierr = ISGetLocalSize(iscol,&ncol);CHKERRQ(ierr); 2184 if (idflag && ncol == mat->cmap->N) allcols = PETSC_TRUE; 2185 else allcols = PETSC_FALSE; 2186 2187 ierr = ISIdentity(isrow,&idflag);CHKERRQ(ierr); 2188 ierr = ISGetLocalSize(isrow,&nrow);CHKERRQ(ierr); 2189 if (idflag && nrow == mat->rmap->N) allrows = PETSC_TRUE; 2190 else allrows = PETSC_FALSE; 2191 2192 if (call == MAT_REUSE_MATRIX) { 2193 ierr = PetscObjectQuery((PetscObject)*newmat,"SubMatrix",(PetscObject*)&Mreuse);CHKERRQ(ierr); 2194 if (!Mreuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Submatrix passed in was not used before, cannot reuse"); 2195 ierr = MatGetSubMatrices_MPIBAIJ_local(mat,1,&isrow_new,&iscol_new,MAT_REUSE_MATRIX,&allrows,&allcols,&Mreuse);CHKERRQ(ierr); 2196 } else { 2197 ierr = MatGetSubMatrices_MPIBAIJ_local(mat,1,&isrow_new,&iscol_new,MAT_INITIAL_MATRIX,&allrows,&allcols,&Mreuse);CHKERRQ(ierr); 2198 } 2199 ierr = ISDestroy(&isrow_new);CHKERRQ(ierr); 2200 ierr = ISDestroy(&iscol_new);CHKERRQ(ierr); 2201 /* 2202 m - number of local rows 2203 n - number of columns (same on all processors) 2204 rstart - first row in new global matrix generated 2205 */ 2206 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 2207 ierr = MatGetSize(Mreuse,&m,&n);CHKERRQ(ierr); 2208 m = m/bs; 2209 n = n/bs; 2210 2211 if (call == MAT_INITIAL_MATRIX) { 2212 aij = (Mat_SeqBAIJ*)(Mreuse)->data; 2213 ii = aij->i; 2214 jj = aij->j; 2215 2216 /* 2217 Determine the number of non-zeros in the diagonal and off-diagonal 2218 portions of the matrix in order to do correct preallocation 2219 */ 2220 2221 /* first get start and end of "diagonal" columns */ 2222 if (csize == PETSC_DECIDE) { 2223 ierr = ISGetSize(isrow,&mglobal);CHKERRQ(ierr); 2224 if (mglobal == n*bs) { /* square matrix */ 2225 nlocal = m; 2226 } else { 2227 nlocal = n/size + ((n % size) > rank); 2228 } 2229 } else { 2230 nlocal = csize/bs; 2231 } 2232 ierr = MPI_Scan(&nlocal,&rend,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); 2233 rstart = rend - nlocal; 2234 if (rank == size - 1 && rend != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local column sizes %D do not add up to total number of columns %D",rend,n); 2235 2236 /* next, compute all the lengths */ 2237 ierr = PetscMalloc2(m+1,&dlens,m+1,&olens);CHKERRQ(ierr); 2238 for (i=0; i<m; i++) { 2239 jend = ii[i+1] - ii[i]; 2240 olen = 0; 2241 dlen = 0; 2242 for (j=0; j<jend; j++) { 2243 if (*jj < rstart || *jj >= rend) olen++; 2244 else dlen++; 2245 jj++; 2246 } 2247 olens[i] = olen; 2248 dlens[i] = dlen; 2249 } 2250 ierr = MatCreate(comm,&M);CHKERRQ(ierr); 2251 ierr = MatSetSizes(M,bs*m,bs*nlocal,PETSC_DECIDE,bs*n);CHKERRQ(ierr); 2252 ierr = MatSetType(M,((PetscObject)mat)->type_name);CHKERRQ(ierr); 2253 ierr = MatMPIBAIJSetPreallocation(M,bs,0,dlens,0,olens);CHKERRQ(ierr); 2254 ierr = MatMPISBAIJSetPreallocation(M,bs,0,dlens,0,olens);CHKERRQ(ierr); 2255 ierr = PetscFree2(dlens,olens);CHKERRQ(ierr); 2256 } else { 2257 PetscInt ml,nl; 2258 2259 M = *newmat; 2260 ierr = MatGetLocalSize(M,&ml,&nl);CHKERRQ(ierr); 2261 if (ml != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Previous matrix must be same size/layout as request"); 2262 ierr = MatZeroEntries(M);CHKERRQ(ierr); 2263 /* 2264 The next two lines are needed so we may call MatSetValues_MPIAIJ() below directly, 2265 rather than the slower MatSetValues(). 2266 */ 2267 M->was_assembled = PETSC_TRUE; 2268 M->assembled = PETSC_FALSE; 2269 } 2270 ierr = MatSetOption(M,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 2271 ierr = MatGetOwnershipRange(M,&rstart,&rend);CHKERRQ(ierr); 2272 aij = (Mat_SeqBAIJ*)(Mreuse)->data; 2273 ii = aij->i; 2274 jj = aij->j; 2275 aa = aij->a; 2276 for (i=0; i<m; i++) { 2277 row = rstart/bs + i; 2278 nz = ii[i+1] - ii[i]; 2279 cwork = jj; jj += nz; 2280 vwork = aa; aa += nz*bs*bs; 2281 ierr = MatSetValuesBlocked_MPIBAIJ(M,1,&row,nz,cwork,vwork,INSERT_VALUES);CHKERRQ(ierr); 2282 } 2283 2284 ierr = MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2285 ierr = MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2286 *newmat = M; 2287 2288 /* save submatrix used in processor for next request */ 2289 if (call == MAT_INITIAL_MATRIX) { 2290 ierr = PetscObjectCompose((PetscObject)M,"SubMatrix",(PetscObject)Mreuse);CHKERRQ(ierr); 2291 ierr = PetscObjectDereference((PetscObject)Mreuse);CHKERRQ(ierr); 2292 } 2293 PetscFunctionReturn(0); 2294 } 2295 2296 #undef __FUNCT__ 2297 #define __FUNCT__ "MatPermute_MPIBAIJ" 2298 PetscErrorCode MatPermute_MPIBAIJ(Mat A,IS rowp,IS colp,Mat *B) 2299 { 2300 MPI_Comm comm,pcomm; 2301 PetscInt clocal_size,nrows; 2302 const PetscInt *rows; 2303 PetscMPIInt size; 2304 IS crowp,lcolp; 2305 PetscErrorCode ierr; 2306 2307 PetscFunctionBegin; 2308 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 2309 /* make a collective version of 'rowp' */ 2310 ierr = PetscObjectGetComm((PetscObject)rowp,&pcomm);CHKERRQ(ierr); 2311 if (pcomm==comm) { 2312 crowp = rowp; 2313 } else { 2314 ierr = ISGetSize(rowp,&nrows);CHKERRQ(ierr); 2315 ierr = ISGetIndices(rowp,&rows);CHKERRQ(ierr); 2316 ierr = ISCreateGeneral(comm,nrows,rows,PETSC_COPY_VALUES,&crowp);CHKERRQ(ierr); 2317 ierr = ISRestoreIndices(rowp,&rows);CHKERRQ(ierr); 2318 } 2319 ierr = ISSetPermutation(crowp);CHKERRQ(ierr); 2320 /* make a local version of 'colp' */ 2321 ierr = PetscObjectGetComm((PetscObject)colp,&pcomm);CHKERRQ(ierr); 2322 ierr = MPI_Comm_size(pcomm,&size);CHKERRQ(ierr); 2323 if (size==1) { 2324 lcolp = colp; 2325 } else { 2326 ierr = ISAllGather(colp,&lcolp);CHKERRQ(ierr); 2327 } 2328 ierr = ISSetPermutation(lcolp);CHKERRQ(ierr); 2329 /* now we just get the submatrix */ 2330 ierr = MatGetLocalSize(A,NULL,&clocal_size);CHKERRQ(ierr); 2331 ierr = MatGetSubMatrix_MPIBAIJ_Private(A,crowp,lcolp,clocal_size,MAT_INITIAL_MATRIX,B);CHKERRQ(ierr); 2332 /* clean up */ 2333 if (pcomm!=comm) { 2334 ierr = ISDestroy(&crowp);CHKERRQ(ierr); 2335 } 2336 if (size>1) { 2337 ierr = ISDestroy(&lcolp);CHKERRQ(ierr); 2338 } 2339 PetscFunctionReturn(0); 2340 } 2341 2342 #undef __FUNCT__ 2343 #define __FUNCT__ "MatGetGhosts_MPIBAIJ" 2344 PetscErrorCode MatGetGhosts_MPIBAIJ(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[]) 2345 { 2346 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*) mat->data; 2347 Mat_SeqBAIJ *B = (Mat_SeqBAIJ*)baij->B->data; 2348 2349 PetscFunctionBegin; 2350 if (nghosts) *nghosts = B->nbs; 2351 if (ghosts) *ghosts = baij->garray; 2352 PetscFunctionReturn(0); 2353 } 2354 2355 #undef __FUNCT__ 2356 #define __FUNCT__ "MatGetSeqNonzeroStructure_MPIBAIJ" 2357 PetscErrorCode MatGetSeqNonzeroStructure_MPIBAIJ(Mat A,Mat *newmat) 2358 { 2359 Mat B; 2360 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 2361 Mat_SeqBAIJ *ad = (Mat_SeqBAIJ*)a->A->data,*bd = (Mat_SeqBAIJ*)a->B->data; 2362 Mat_SeqAIJ *b; 2363 PetscErrorCode ierr; 2364 PetscMPIInt size,rank,*recvcounts = 0,*displs = 0; 2365 PetscInt sendcount,i,*rstarts = A->rmap->range,n,cnt,j,bs = A->rmap->bs; 2366 PetscInt m,*garray = a->garray,*lens,*jsendbuf,*a_jsendbuf,*b_jsendbuf; 2367 2368 PetscFunctionBegin; 2369 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 2370 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)A),&rank);CHKERRQ(ierr); 2371 2372 /* ---------------------------------------------------------------- 2373 Tell every processor the number of nonzeros per row 2374 */ 2375 ierr = PetscMalloc1(A->rmap->N/bs,&lens);CHKERRQ(ierr); 2376 for (i=A->rmap->rstart/bs; i<A->rmap->rend/bs; i++) { 2377 lens[i] = ad->i[i-A->rmap->rstart/bs+1] - ad->i[i-A->rmap->rstart/bs] + bd->i[i-A->rmap->rstart/bs+1] - bd->i[i-A->rmap->rstart/bs]; 2378 } 2379 ierr = PetscMalloc1(2*size,&recvcounts);CHKERRQ(ierr); 2380 displs = recvcounts + size; 2381 for (i=0; i<size; i++) { 2382 recvcounts[i] = A->rmap->range[i+1]/bs - A->rmap->range[i]/bs; 2383 displs[i] = A->rmap->range[i]/bs; 2384 } 2385 #if defined(PETSC_HAVE_MPI_IN_PLACE) 2386 ierr = MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,lens,recvcounts,displs,MPIU_INT,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2387 #else 2388 sendcount = A->rmap->rend/bs - A->rmap->rstart/bs; 2389 ierr = MPI_Allgatherv(lens+A->rmap->rstart/bs,sendcount,MPIU_INT,lens,recvcounts,displs,MPIU_INT,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2390 #endif 2391 /* --------------------------------------------------------------- 2392 Create the sequential matrix of the same type as the local block diagonal 2393 */ 2394 ierr = MatCreate(PETSC_COMM_SELF,&B);CHKERRQ(ierr); 2395 ierr = MatSetSizes(B,A->rmap->N/bs,A->cmap->N/bs,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2396 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2397 ierr = MatSeqAIJSetPreallocation(B,0,lens);CHKERRQ(ierr); 2398 b = (Mat_SeqAIJ*)B->data; 2399 2400 /*-------------------------------------------------------------------- 2401 Copy my part of matrix column indices over 2402 */ 2403 sendcount = ad->nz + bd->nz; 2404 jsendbuf = b->j + b->i[rstarts[rank]/bs]; 2405 a_jsendbuf = ad->j; 2406 b_jsendbuf = bd->j; 2407 n = A->rmap->rend/bs - A->rmap->rstart/bs; 2408 cnt = 0; 2409 for (i=0; i<n; i++) { 2410 2411 /* put in lower diagonal portion */ 2412 m = bd->i[i+1] - bd->i[i]; 2413 while (m > 0) { 2414 /* is it above diagonal (in bd (compressed) numbering) */ 2415 if (garray[*b_jsendbuf] > A->rmap->rstart/bs + i) break; 2416 jsendbuf[cnt++] = garray[*b_jsendbuf++]; 2417 m--; 2418 } 2419 2420 /* put in diagonal portion */ 2421 for (j=ad->i[i]; j<ad->i[i+1]; j++) { 2422 jsendbuf[cnt++] = A->rmap->rstart/bs + *a_jsendbuf++; 2423 } 2424 2425 /* put in upper diagonal portion */ 2426 while (m-- > 0) { 2427 jsendbuf[cnt++] = garray[*b_jsendbuf++]; 2428 } 2429 } 2430 if (cnt != sendcount) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupted PETSc matrix: nz given %D actual nz %D",sendcount,cnt); 2431 2432 /*-------------------------------------------------------------------- 2433 Gather all column indices to all processors 2434 */ 2435 for (i=0; i<size; i++) { 2436 recvcounts[i] = 0; 2437 for (j=A->rmap->range[i]/bs; j<A->rmap->range[i+1]/bs; j++) { 2438 recvcounts[i] += lens[j]; 2439 } 2440 } 2441 displs[0] = 0; 2442 for (i=1; i<size; i++) { 2443 displs[i] = displs[i-1] + recvcounts[i-1]; 2444 } 2445 #if defined(PETSC_HAVE_MPI_IN_PLACE) 2446 ierr = MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,b->j,recvcounts,displs,MPIU_INT,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2447 #else 2448 ierr = MPI_Allgatherv(jsendbuf,sendcount,MPIU_INT,b->j,recvcounts,displs,MPIU_INT,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2449 #endif 2450 /*-------------------------------------------------------------------- 2451 Assemble the matrix into useable form (note numerical values not yet set) 2452 */ 2453 /* set the b->ilen (length of each row) values */ 2454 ierr = PetscMemcpy(b->ilen,lens,(A->rmap->N/bs)*sizeof(PetscInt));CHKERRQ(ierr); 2455 /* set the b->i indices */ 2456 b->i[0] = 0; 2457 for (i=1; i<=A->rmap->N/bs; i++) { 2458 b->i[i] = b->i[i-1] + lens[i-1]; 2459 } 2460 ierr = PetscFree(lens);CHKERRQ(ierr); 2461 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2462 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2463 ierr = PetscFree(recvcounts);CHKERRQ(ierr); 2464 2465 if (A->symmetric) { 2466 ierr = MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 2467 } else if (A->hermitian) { 2468 ierr = MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE);CHKERRQ(ierr); 2469 } else if (A->structurally_symmetric) { 2470 ierr = MatSetOption(B,MAT_STRUCTURALLY_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 2471 } 2472 *newmat = B; 2473 PetscFunctionReturn(0); 2474 } 2475 2476 #undef __FUNCT__ 2477 #define __FUNCT__ "MatSOR_MPIBAIJ" 2478 PetscErrorCode MatSOR_MPIBAIJ(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 2479 { 2480 Mat_MPIBAIJ *mat = (Mat_MPIBAIJ*)matin->data; 2481 PetscErrorCode ierr; 2482 Vec bb1 = 0; 2483 2484 PetscFunctionBegin; 2485 if (flag == SOR_APPLY_UPPER) { 2486 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 2487 PetscFunctionReturn(0); 2488 } 2489 2490 if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS) { 2491 ierr = VecDuplicate(bb,&bb1);CHKERRQ(ierr); 2492 } 2493 2494 if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) { 2495 if (flag & SOR_ZERO_INITIAL_GUESS) { 2496 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 2497 its--; 2498 } 2499 2500 while (its--) { 2501 ierr = VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2502 ierr = VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2503 2504 /* update rhs: bb1 = bb - B*x */ 2505 ierr = VecScale(mat->lvec,-1.0);CHKERRQ(ierr); 2506 ierr = (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);CHKERRQ(ierr); 2507 2508 /* local sweep */ 2509 ierr = (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx);CHKERRQ(ierr); 2510 } 2511 } else if (flag & SOR_LOCAL_FORWARD_SWEEP) { 2512 if (flag & SOR_ZERO_INITIAL_GUESS) { 2513 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 2514 its--; 2515 } 2516 while (its--) { 2517 ierr = VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2518 ierr = VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2519 2520 /* update rhs: bb1 = bb - B*x */ 2521 ierr = VecScale(mat->lvec,-1.0);CHKERRQ(ierr); 2522 ierr = (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);CHKERRQ(ierr); 2523 2524 /* local sweep */ 2525 ierr = (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx);CHKERRQ(ierr); 2526 } 2527 } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) { 2528 if (flag & SOR_ZERO_INITIAL_GUESS) { 2529 ierr = (*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx);CHKERRQ(ierr); 2530 its--; 2531 } 2532 while (its--) { 2533 ierr = VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2534 ierr = VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2535 2536 /* update rhs: bb1 = bb - B*x */ 2537 ierr = VecScale(mat->lvec,-1.0);CHKERRQ(ierr); 2538 ierr = (*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1);CHKERRQ(ierr); 2539 2540 /* local sweep */ 2541 ierr = (*mat->A->ops->sor)(mat->A,bb1,omega,SOR_BACKWARD_SWEEP,fshift,lits,1,xx);CHKERRQ(ierr); 2542 } 2543 } else SETERRQ(PetscObjectComm((PetscObject)matin),PETSC_ERR_SUP,"Parallel version of SOR requested not supported"); 2544 2545 ierr = VecDestroy(&bb1);CHKERRQ(ierr); 2546 PetscFunctionReturn(0); 2547 } 2548 2549 #undef __FUNCT__ 2550 #define __FUNCT__ "MatGetColumnNorms_MPIBAIJ" 2551 PetscErrorCode MatGetColumnNorms_MPIBAIJ(Mat A,NormType type,PetscReal *norms) 2552 { 2553 PetscErrorCode ierr; 2554 Mat_MPIBAIJ *aij = (Mat_MPIBAIJ*)A->data; 2555 PetscInt N,i,*garray = aij->garray; 2556 PetscInt ib,jb,bs = A->rmap->bs; 2557 Mat_SeqBAIJ *a_aij = (Mat_SeqBAIJ*) aij->A->data; 2558 MatScalar *a_val = a_aij->a; 2559 Mat_SeqBAIJ *b_aij = (Mat_SeqBAIJ*) aij->B->data; 2560 MatScalar *b_val = b_aij->a; 2561 PetscReal *work; 2562 2563 PetscFunctionBegin; 2564 ierr = MatGetSize(A,NULL,&N);CHKERRQ(ierr); 2565 ierr = PetscCalloc1(N,&work);CHKERRQ(ierr); 2566 if (type == NORM_2) { 2567 for (i=a_aij->i[0]; i<a_aij->i[aij->A->rmap->n/bs]; i++) { 2568 for (jb=0; jb<bs; jb++) { 2569 for (ib=0; ib<bs; ib++) { 2570 work[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscAbsScalar(*a_val * *a_val); 2571 a_val++; 2572 } 2573 } 2574 } 2575 for (i=b_aij->i[0]; i<b_aij->i[aij->B->rmap->n/bs]; i++) { 2576 for (jb=0; jb<bs; jb++) { 2577 for (ib=0; ib<bs; ib++) { 2578 work[garray[b_aij->j[i]] * bs + jb] += PetscAbsScalar(*b_val * *b_val); 2579 b_val++; 2580 } 2581 } 2582 } 2583 } else if (type == NORM_1) { 2584 for (i=a_aij->i[0]; i<a_aij->i[aij->A->rmap->n/bs]; i++) { 2585 for (jb=0; jb<bs; jb++) { 2586 for (ib=0; ib<bs; ib++) { 2587 work[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscAbsScalar(*a_val); 2588 a_val++; 2589 } 2590 } 2591 } 2592 for (i=b_aij->i[0]; i<b_aij->i[aij->B->rmap->n/bs]; i++) { 2593 for (jb=0; jb<bs; jb++) { 2594 for (ib=0; ib<bs; ib++) { 2595 work[garray[b_aij->j[i]] * bs + jb] += PetscAbsScalar(*b_val); 2596 b_val++; 2597 } 2598 } 2599 } 2600 } else if (type == NORM_INFINITY) { 2601 for (i=a_aij->i[0]; i<a_aij->i[aij->A->rmap->n/bs]; i++) { 2602 for (jb=0; jb<bs; jb++) { 2603 for (ib=0; ib<bs; ib++) { 2604 int col = A->cmap->rstart + a_aij->j[i] * bs + jb; 2605 work[col] = PetscMax(PetscAbsScalar(*a_val), work[col]); 2606 a_val++; 2607 } 2608 } 2609 } 2610 for (i=b_aij->i[0]; i<b_aij->i[aij->B->rmap->n/bs]; i++) { 2611 for (jb=0; jb<bs; jb++) { 2612 for (ib=0; ib<bs; ib++) { 2613 int col = garray[b_aij->j[i]] * bs + jb; 2614 work[col] = PetscMax(PetscAbsScalar(*b_val), work[col]); 2615 b_val++; 2616 } 2617 } 2618 } 2619 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2620 if (type == NORM_INFINITY) { 2621 ierr = MPIU_Allreduce(work,norms,N,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2622 } else { 2623 ierr = MPIU_Allreduce(work,norms,N,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2624 } 2625 ierr = PetscFree(work);CHKERRQ(ierr); 2626 if (type == NORM_2) { 2627 for (i=0; i<N; i++) norms[i] = PetscSqrtReal(norms[i]); 2628 } 2629 PetscFunctionReturn(0); 2630 } 2631 2632 #undef __FUNCT__ 2633 #define __FUNCT__ "MatInvertBlockDiagonal_MPIBAIJ" 2634 PetscErrorCode MatInvertBlockDiagonal_MPIBAIJ(Mat A,const PetscScalar **values) 2635 { 2636 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*) A->data; 2637 PetscErrorCode ierr; 2638 2639 PetscFunctionBegin; 2640 ierr = MatInvertBlockDiagonal(a->A,values);CHKERRQ(ierr); 2641 A->factorerrortype = a->A->factorerrortype; 2642 A->factorerror_zeropivot_value = a->A->factorerror_zeropivot_value; 2643 A->factorerror_zeropivot_row = a->A->factorerror_zeropivot_row; 2644 PetscFunctionReturn(0); 2645 } 2646 2647 #undef __FUNCT__ 2648 #define __FUNCT__ "MatShift_MPIBAIJ" 2649 PetscErrorCode MatShift_MPIBAIJ(Mat Y,PetscScalar a) 2650 { 2651 PetscErrorCode ierr; 2652 Mat_MPIBAIJ *maij = (Mat_MPIBAIJ*)Y->data; 2653 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ*)maij->A->data; 2654 2655 PetscFunctionBegin; 2656 if (!Y->preallocated) { 2657 ierr = MatMPIBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL,0,NULL);CHKERRQ(ierr); 2658 } else if (!aij->nz) { 2659 PetscInt nonew = aij->nonew; 2660 ierr = MatSeqBAIJSetPreallocation(maij->A,Y->rmap->bs,1,NULL);CHKERRQ(ierr); 2661 aij->nonew = nonew; 2662 } 2663 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 2664 PetscFunctionReturn(0); 2665 } 2666 2667 #undef __FUNCT__ 2668 #define __FUNCT__ "MatMissingDiagonal_MPIBAIJ" 2669 PetscErrorCode MatMissingDiagonal_MPIBAIJ(Mat A,PetscBool *missing,PetscInt *d) 2670 { 2671 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 2672 PetscErrorCode ierr; 2673 2674 PetscFunctionBegin; 2675 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only works for square matrices"); 2676 ierr = MatMissingDiagonal(a->A,missing,d);CHKERRQ(ierr); 2677 if (d) { 2678 PetscInt rstart; 2679 ierr = MatGetOwnershipRange(A,&rstart,NULL);CHKERRQ(ierr); 2680 *d += rstart/A->rmap->bs; 2681 2682 } 2683 PetscFunctionReturn(0); 2684 } 2685 2686 #undef __FUNCT__ 2687 #define __FUNCT__ "MatGetDiagonalBlock_MPIBAIJ" 2688 PetscErrorCode MatGetDiagonalBlock_MPIBAIJ(Mat A,Mat *a) 2689 { 2690 PetscFunctionBegin; 2691 *a = ((Mat_MPIBAIJ*)A->data)->A; 2692 PetscFunctionReturn(0); 2693 } 2694 2695 /* -------------------------------------------------------------------*/ 2696 static struct _MatOps MatOps_Values = {MatSetValues_MPIBAIJ, 2697 MatGetRow_MPIBAIJ, 2698 MatRestoreRow_MPIBAIJ, 2699 MatMult_MPIBAIJ, 2700 /* 4*/ MatMultAdd_MPIBAIJ, 2701 MatMultTranspose_MPIBAIJ, 2702 MatMultTransposeAdd_MPIBAIJ, 2703 0, 2704 0, 2705 0, 2706 /*10*/ 0, 2707 0, 2708 0, 2709 MatSOR_MPIBAIJ, 2710 MatTranspose_MPIBAIJ, 2711 /*15*/ MatGetInfo_MPIBAIJ, 2712 MatEqual_MPIBAIJ, 2713 MatGetDiagonal_MPIBAIJ, 2714 MatDiagonalScale_MPIBAIJ, 2715 MatNorm_MPIBAIJ, 2716 /*20*/ MatAssemblyBegin_MPIBAIJ, 2717 MatAssemblyEnd_MPIBAIJ, 2718 MatSetOption_MPIBAIJ, 2719 MatZeroEntries_MPIBAIJ, 2720 /*24*/ MatZeroRows_MPIBAIJ, 2721 0, 2722 0, 2723 0, 2724 0, 2725 /*29*/ MatSetUp_MPIBAIJ, 2726 0, 2727 0, 2728 MatGetDiagonalBlock_MPIBAIJ, 2729 0, 2730 /*34*/ MatDuplicate_MPIBAIJ, 2731 0, 2732 0, 2733 0, 2734 0, 2735 /*39*/ MatAXPY_MPIBAIJ, 2736 MatGetSubMatrices_MPIBAIJ, 2737 MatIncreaseOverlap_MPIBAIJ, 2738 MatGetValues_MPIBAIJ, 2739 MatCopy_MPIBAIJ, 2740 /*44*/ 0, 2741 MatScale_MPIBAIJ, 2742 MatShift_MPIBAIJ, 2743 0, 2744 MatZeroRowsColumns_MPIBAIJ, 2745 /*49*/ 0, 2746 0, 2747 0, 2748 0, 2749 0, 2750 /*54*/ MatFDColoringCreate_MPIXAIJ, 2751 0, 2752 MatSetUnfactored_MPIBAIJ, 2753 MatPermute_MPIBAIJ, 2754 MatSetValuesBlocked_MPIBAIJ, 2755 /*59*/ MatGetSubMatrix_MPIBAIJ, 2756 MatDestroy_MPIBAIJ, 2757 MatView_MPIBAIJ, 2758 0, 2759 0, 2760 /*64*/ 0, 2761 0, 2762 0, 2763 0, 2764 0, 2765 /*69*/ MatGetRowMaxAbs_MPIBAIJ, 2766 0, 2767 0, 2768 0, 2769 0, 2770 /*74*/ 0, 2771 MatFDColoringApply_BAIJ, 2772 0, 2773 0, 2774 0, 2775 /*79*/ 0, 2776 0, 2777 0, 2778 0, 2779 MatLoad_MPIBAIJ, 2780 /*84*/ 0, 2781 0, 2782 0, 2783 0, 2784 0, 2785 /*89*/ 0, 2786 0, 2787 0, 2788 0, 2789 0, 2790 /*94*/ 0, 2791 0, 2792 0, 2793 0, 2794 0, 2795 /*99*/ 0, 2796 0, 2797 0, 2798 0, 2799 0, 2800 /*104*/0, 2801 MatRealPart_MPIBAIJ, 2802 MatImaginaryPart_MPIBAIJ, 2803 0, 2804 0, 2805 /*109*/0, 2806 0, 2807 0, 2808 0, 2809 MatMissingDiagonal_MPIBAIJ, 2810 /*114*/MatGetSeqNonzeroStructure_MPIBAIJ, 2811 0, 2812 MatGetGhosts_MPIBAIJ, 2813 0, 2814 0, 2815 /*119*/0, 2816 0, 2817 0, 2818 0, 2819 MatGetMultiProcBlock_MPIBAIJ, 2820 /*124*/0, 2821 MatGetColumnNorms_MPIBAIJ, 2822 MatInvertBlockDiagonal_MPIBAIJ, 2823 0, 2824 0, 2825 /*129*/ 0, 2826 0, 2827 0, 2828 0, 2829 0, 2830 /*134*/ 0, 2831 0, 2832 0, 2833 0, 2834 0, 2835 /*139*/ 0, 2836 0, 2837 0, 2838 MatFDColoringSetUp_MPIXAIJ, 2839 0, 2840 /*144*/MatCreateMPIMatConcatenateSeqMat_MPIBAIJ 2841 }; 2842 2843 2844 PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPISBAIJ(Mat, MatType,MatReuse,Mat*); 2845 2846 #undef __FUNCT__ 2847 #define __FUNCT__ "MatMPIBAIJSetPreallocationCSR_MPIBAIJ" 2848 PetscErrorCode MatMPIBAIJSetPreallocationCSR_MPIBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[],const PetscScalar V[]) 2849 { 2850 PetscInt m,rstart,cstart,cend; 2851 PetscInt i,j,d,nz,nz_max=0,*d_nnz=0,*o_nnz=0; 2852 const PetscInt *JJ =0; 2853 PetscScalar *values=0; 2854 PetscBool roworiented = ((Mat_MPIBAIJ*)B->data)->roworiented; 2855 PetscErrorCode ierr; 2856 2857 PetscFunctionBegin; 2858 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 2859 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 2860 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2861 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2862 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 2863 m = B->rmap->n/bs; 2864 rstart = B->rmap->rstart/bs; 2865 cstart = B->cmap->rstart/bs; 2866 cend = B->cmap->rend/bs; 2867 2868 if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]); 2869 ierr = PetscMalloc2(m,&d_nnz,m,&o_nnz);CHKERRQ(ierr); 2870 for (i=0; i<m; i++) { 2871 nz = ii[i+1] - ii[i]; 2872 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local row %D has a negative number of columns %D",i,nz); 2873 nz_max = PetscMax(nz_max,nz); 2874 JJ = jj + ii[i]; 2875 for (j=0; j<nz; j++) { 2876 if (*JJ >= cstart) break; 2877 JJ++; 2878 } 2879 d = 0; 2880 for (; j<nz; j++) { 2881 if (*JJ++ >= cend) break; 2882 d++; 2883 } 2884 d_nnz[i] = d; 2885 o_nnz[i] = nz - d; 2886 } 2887 ierr = MatMPIBAIJSetPreallocation(B,bs,0,d_nnz,0,o_nnz);CHKERRQ(ierr); 2888 ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr); 2889 2890 values = (PetscScalar*)V; 2891 if (!values) { 2892 ierr = PetscMalloc1(bs*bs*nz_max,&values);CHKERRQ(ierr); 2893 ierr = PetscMemzero(values,bs*bs*nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 2894 } 2895 for (i=0; i<m; i++) { 2896 PetscInt row = i + rstart; 2897 PetscInt ncols = ii[i+1] - ii[i]; 2898 const PetscInt *icols = jj + ii[i]; 2899 if (!roworiented) { /* block ordering matches the non-nested layout of MatSetValues so we can insert entire rows */ 2900 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 2901 ierr = MatSetValuesBlocked_MPIBAIJ(B,1,&row,ncols,icols,svals,INSERT_VALUES);CHKERRQ(ierr); 2902 } else { /* block ordering does not match so we can only insert one block at a time. */ 2903 PetscInt j; 2904 for (j=0; j<ncols; j++) { 2905 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 2906 ierr = MatSetValuesBlocked_MPIBAIJ(B,1,&row,1,&icols[j],svals,INSERT_VALUES);CHKERRQ(ierr); 2907 } 2908 } 2909 } 2910 2911 if (!V) { ierr = PetscFree(values);CHKERRQ(ierr); } 2912 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2913 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2914 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2915 PetscFunctionReturn(0); 2916 } 2917 2918 #undef __FUNCT__ 2919 #define __FUNCT__ "MatMPIBAIJSetPreallocationCSR" 2920 /*@C 2921 MatMPIBAIJSetPreallocationCSR - Allocates memory for a sparse parallel matrix in BAIJ format 2922 (the default parallel PETSc format). 2923 2924 Collective on MPI_Comm 2925 2926 Input Parameters: 2927 + B - the matrix 2928 . bs - the block size 2929 . i - the indices into j for the start of each local row (starts with zero) 2930 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2931 - v - optional values in the matrix 2932 2933 Level: developer 2934 2935 Notes: The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2936 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2937 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2938 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2939 block column and the second index is over columns within a block. 2940 2941 .keywords: matrix, aij, compressed row, sparse, parallel 2942 2943 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatMPIBAIJSetPreallocation(), MatCreateAIJ(), MPIAIJ, MatCreateMPIBAIJWithArrays(), MPIBAIJ 2944 @*/ 2945 PetscErrorCode MatMPIBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2946 { 2947 PetscErrorCode ierr; 2948 2949 PetscFunctionBegin; 2950 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2951 PetscValidType(B,1); 2952 PetscValidLogicalCollectiveInt(B,bs,2); 2953 ierr = PetscTryMethod(B,"MatMPIBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));CHKERRQ(ierr); 2954 PetscFunctionReturn(0); 2955 } 2956 2957 #undef __FUNCT__ 2958 #define __FUNCT__ "MatMPIBAIJSetPreallocation_MPIBAIJ" 2959 PetscErrorCode MatMPIBAIJSetPreallocation_MPIBAIJ(Mat B,PetscInt bs,PetscInt d_nz,const PetscInt *d_nnz,PetscInt o_nz,const PetscInt *o_nnz) 2960 { 2961 Mat_MPIBAIJ *b; 2962 PetscErrorCode ierr; 2963 PetscInt i; 2964 2965 PetscFunctionBegin; 2966 ierr = MatSetBlockSize(B,PetscAbs(bs));CHKERRQ(ierr); 2967 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2968 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2969 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 2970 2971 if (d_nnz) { 2972 for (i=0; i<B->rmap->n/bs; i++) { 2973 if (d_nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"d_nnz cannot be less than -1: local row %D value %D",i,d_nnz[i]); 2974 } 2975 } 2976 if (o_nnz) { 2977 for (i=0; i<B->rmap->n/bs; i++) { 2978 if (o_nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"o_nnz cannot be less than -1: local row %D value %D",i,o_nnz[i]); 2979 } 2980 } 2981 2982 b = (Mat_MPIBAIJ*)B->data; 2983 b->bs2 = bs*bs; 2984 b->mbs = B->rmap->n/bs; 2985 b->nbs = B->cmap->n/bs; 2986 b->Mbs = B->rmap->N/bs; 2987 b->Nbs = B->cmap->N/bs; 2988 2989 for (i=0; i<=b->size; i++) { 2990 b->rangebs[i] = B->rmap->range[i]/bs; 2991 } 2992 b->rstartbs = B->rmap->rstart/bs; 2993 b->rendbs = B->rmap->rend/bs; 2994 b->cstartbs = B->cmap->rstart/bs; 2995 b->cendbs = B->cmap->rend/bs; 2996 2997 #if defined(PETSC_USE_CTABLE) 2998 ierr = PetscTableDestroy(&b->colmap);CHKERRQ(ierr); 2999 #else 3000 ierr = PetscFree(b->colmap);CHKERRQ(ierr); 3001 #endif 3002 ierr = PetscFree(b->garray);CHKERRQ(ierr); 3003 ierr = VecDestroy(&b->lvec);CHKERRQ(ierr); 3004 ierr = VecScatterDestroy(&b->Mvctx);CHKERRQ(ierr); 3005 3006 /* Because the B will have been resized we simply destroy it and create a new one each time */ 3007 ierr = MatDestroy(&b->B);CHKERRQ(ierr); 3008 ierr = MatCreate(PETSC_COMM_SELF,&b->B);CHKERRQ(ierr); 3009 ierr = MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N);CHKERRQ(ierr); 3010 ierr = MatSetType(b->B,MATSEQBAIJ);CHKERRQ(ierr); 3011 ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)b->B);CHKERRQ(ierr); 3012 3013 if (!B->preallocated) { 3014 ierr = MatCreate(PETSC_COMM_SELF,&b->A);CHKERRQ(ierr); 3015 ierr = MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n);CHKERRQ(ierr); 3016 ierr = MatSetType(b->A,MATSEQBAIJ);CHKERRQ(ierr); 3017 ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)b->A);CHKERRQ(ierr); 3018 ierr = MatStashCreate_Private(PetscObjectComm((PetscObject)B),bs,&B->bstash);CHKERRQ(ierr); 3019 } 3020 3021 ierr = MatSeqBAIJSetPreallocation(b->A,bs,d_nz,d_nnz);CHKERRQ(ierr); 3022 ierr = MatSeqBAIJSetPreallocation(b->B,bs,o_nz,o_nnz);CHKERRQ(ierr); 3023 B->preallocated = PETSC_TRUE; 3024 B->was_assembled = PETSC_FALSE; 3025 B->assembled = PETSC_FALSE; 3026 PetscFunctionReturn(0); 3027 } 3028 3029 extern PetscErrorCode MatDiagonalScaleLocal_MPIBAIJ(Mat,Vec); 3030 extern PetscErrorCode MatSetHashTableFactor_MPIBAIJ(Mat,PetscReal); 3031 3032 #undef __FUNCT__ 3033 #define __FUNCT__ "MatConvert_MPIBAIJ_MPIAdj" 3034 PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPIAdj(Mat B, MatType newtype,MatReuse reuse,Mat *adj) 3035 { 3036 Mat_MPIBAIJ *b = (Mat_MPIBAIJ*)B->data; 3037 PetscErrorCode ierr; 3038 Mat_SeqBAIJ *d = (Mat_SeqBAIJ*) b->A->data,*o = (Mat_SeqBAIJ*) b->B->data; 3039 PetscInt M = B->rmap->n/B->rmap->bs,i,*ii,*jj,cnt,j,k,rstart = B->rmap->rstart/B->rmap->bs; 3040 const PetscInt *id = d->i, *jd = d->j, *io = o->i, *jo = o->j, *garray = b->garray; 3041 3042 PetscFunctionBegin; 3043 ierr = PetscMalloc1(M+1,&ii);CHKERRQ(ierr); 3044 ii[0] = 0; 3045 for (i=0; i<M; i++) { 3046 if ((id[i+1] - id[i]) < 0) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Indices wrong %D %D %D",i,id[i],id[i+1]); 3047 if ((io[i+1] - io[i]) < 0) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Indices wrong %D %D %D",i,io[i],io[i+1]); 3048 ii[i+1] = ii[i] + id[i+1] - id[i] + io[i+1] - io[i]; 3049 /* remove one from count of matrix has diagonal */ 3050 for (j=id[i]; j<id[i+1]; j++) { 3051 if (jd[j] == i) {ii[i+1]--;break;} 3052 } 3053 } 3054 ierr = PetscMalloc1(ii[M],&jj);CHKERRQ(ierr); 3055 cnt = 0; 3056 for (i=0; i<M; i++) { 3057 for (j=io[i]; j<io[i+1]; j++) { 3058 if (garray[jo[j]] > rstart) break; 3059 jj[cnt++] = garray[jo[j]]; 3060 } 3061 for (k=id[i]; k<id[i+1]; k++) { 3062 if (jd[k] != i) { 3063 jj[cnt++] = rstart + jd[k]; 3064 } 3065 } 3066 for (; j<io[i+1]; j++) { 3067 jj[cnt++] = garray[jo[j]]; 3068 } 3069 } 3070 ierr = MatCreateMPIAdj(PetscObjectComm((PetscObject)B),M,B->cmap->N/B->rmap->bs,ii,jj,NULL,adj);CHKERRQ(ierr); 3071 PetscFunctionReturn(0); 3072 } 3073 3074 #include <../src/mat/impls/aij/mpi/mpiaij.h> 3075 3076 PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat,MatType,MatReuse,Mat*); 3077 3078 #undef __FUNCT__ 3079 #define __FUNCT__ "MatConvert_MPIBAIJ_MPIAIJ" 3080 PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 3081 { 3082 PetscErrorCode ierr; 3083 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 3084 Mat B; 3085 Mat_MPIAIJ *b; 3086 3087 PetscFunctionBegin; 3088 if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 3089 3090 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 3091 ierr = MatSetType(B,MATMPIAIJ);CHKERRQ(ierr); 3092 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 3093 ierr = MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 3094 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 3095 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 3096 b = (Mat_MPIAIJ*) B->data; 3097 3098 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 3099 ierr = MatDestroy(&b->B);CHKERRQ(ierr); 3100 ierr = MatDisAssemble_MPIBAIJ(A);CHKERRQ(ierr); 3101 ierr = MatConvert_SeqBAIJ_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A);CHKERRQ(ierr); 3102 ierr = MatConvert_SeqBAIJ_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B);CHKERRQ(ierr); 3103 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3104 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3105 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3106 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3107 if (reuse == MAT_INPLACE_MATRIX) { 3108 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 3109 } else { 3110 *newmat = B; 3111 } 3112 PetscFunctionReturn(0); 3113 } 3114 3115 /*MC 3116 MATMPIBAIJ - MATMPIBAIJ = "mpibaij" - A matrix type to be used for distributed block sparse matrices. 3117 3118 Options Database Keys: 3119 + -mat_type mpibaij - sets the matrix type to "mpibaij" during a call to MatSetFromOptions() 3120 . -mat_block_size <bs> - set the blocksize used to store the matrix 3121 - -mat_use_hash_table <fact> 3122 3123 Level: beginner 3124 3125 .seealso: MatCreateMPIBAIJ 3126 M*/ 3127 3128 PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPIBSTRM(Mat,MatType,MatReuse,Mat*); 3129 3130 #undef __FUNCT__ 3131 #define __FUNCT__ "MatCreate_MPIBAIJ" 3132 PETSC_EXTERN PetscErrorCode MatCreate_MPIBAIJ(Mat B) 3133 { 3134 Mat_MPIBAIJ *b; 3135 PetscErrorCode ierr; 3136 PetscBool flg = PETSC_FALSE; 3137 3138 PetscFunctionBegin; 3139 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 3140 B->data = (void*)b; 3141 3142 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3143 B->assembled = PETSC_FALSE; 3144 3145 B->insertmode = NOT_SET_VALUES; 3146 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank);CHKERRQ(ierr); 3147 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&b->size);CHKERRQ(ierr); 3148 3149 /* build local table of row and column ownerships */ 3150 ierr = PetscMalloc1(b->size+1,&b->rangebs);CHKERRQ(ierr); 3151 3152 /* build cache for off array entries formed */ 3153 ierr = MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);CHKERRQ(ierr); 3154 3155 b->donotstash = PETSC_FALSE; 3156 b->colmap = NULL; 3157 b->garray = NULL; 3158 b->roworiented = PETSC_TRUE; 3159 3160 /* stuff used in block assembly */ 3161 b->barray = 0; 3162 3163 /* stuff used for matrix vector multiply */ 3164 b->lvec = 0; 3165 b->Mvctx = 0; 3166 3167 /* stuff for MatGetRow() */ 3168 b->rowindices = 0; 3169 b->rowvalues = 0; 3170 b->getrowactive = PETSC_FALSE; 3171 3172 /* hash table stuff */ 3173 b->ht = 0; 3174 b->hd = 0; 3175 b->ht_size = 0; 3176 b->ht_flag = PETSC_FALSE; 3177 b->ht_fact = 0; 3178 b->ht_total_ct = 0; 3179 b->ht_insert_ct = 0; 3180 3181 /* stuff for MatGetSubMatrices_MPIBAIJ_local() */ 3182 b->ijonly = PETSC_FALSE; 3183 3184 3185 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpibaij_mpiadj_C",MatConvert_MPIBAIJ_MPIAdj);CHKERRQ(ierr); 3186 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpibaij_mpiaij_C",MatConvert_MPIBAIJ_MPIAIJ);CHKERRQ(ierr); 3187 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpibaij_mpisbaij_C",MatConvert_MPIBAIJ_MPISBAIJ);CHKERRQ(ierr); 3188 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPIBAIJ);CHKERRQ(ierr); 3189 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPIBAIJ);CHKERRQ(ierr); 3190 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPIBAIJSetPreallocation_C",MatMPIBAIJSetPreallocation_MPIBAIJ);CHKERRQ(ierr); 3191 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPIBAIJSetPreallocationCSR_C",MatMPIBAIJSetPreallocationCSR_MPIBAIJ);CHKERRQ(ierr); 3192 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPIBAIJ);CHKERRQ(ierr); 3193 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSetHashTableFactor_C",MatSetHashTableFactor_MPIBAIJ);CHKERRQ(ierr); 3194 ierr = PetscObjectChangeTypeName((PetscObject)B,MATMPIBAIJ);CHKERRQ(ierr); 3195 3196 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),NULL,"Options for loading MPIBAIJ matrix 1","Mat");CHKERRQ(ierr); 3197 ierr = PetscOptionsName("-mat_use_hash_table","Use hash table to save time in constructing matrix","MatSetOption",&flg);CHKERRQ(ierr); 3198 if (flg) { 3199 PetscReal fact = 1.39; 3200 ierr = MatSetOption(B,MAT_USE_HASH_TABLE,PETSC_TRUE);CHKERRQ(ierr); 3201 ierr = PetscOptionsReal("-mat_use_hash_table","Use hash table factor","MatMPIBAIJSetHashTableFactor",fact,&fact,NULL);CHKERRQ(ierr); 3202 if (fact <= 1.0) fact = 1.39; 3203 ierr = MatMPIBAIJSetHashTableFactor(B,fact);CHKERRQ(ierr); 3204 ierr = PetscInfo1(B,"Hash table Factor used %5.2f\n",fact);CHKERRQ(ierr); 3205 } 3206 ierr = PetscOptionsEnd();CHKERRQ(ierr); 3207 PetscFunctionReturn(0); 3208 } 3209 3210 /*MC 3211 MATBAIJ - MATBAIJ = "baij" - A matrix type to be used for block sparse matrices. 3212 3213 This matrix type is identical to MATSEQBAIJ when constructed with a single process communicator, 3214 and MATMPIBAIJ otherwise. 3215 3216 Options Database Keys: 3217 . -mat_type baij - sets the matrix type to "baij" during a call to MatSetFromOptions() 3218 3219 Level: beginner 3220 3221 .seealso: MatCreateBAIJ(),MATSEQBAIJ,MATMPIBAIJ, MatMPIBAIJSetPreallocation(), MatMPIBAIJSetPreallocationCSR() 3222 M*/ 3223 3224 #undef __FUNCT__ 3225 #define __FUNCT__ "MatMPIBAIJSetPreallocation" 3226 /*@C 3227 MatMPIBAIJSetPreallocation - Allocates memory for a sparse parallel matrix in block AIJ format 3228 (block compressed row). For good matrix assembly performance 3229 the user should preallocate the matrix storage by setting the parameters 3230 d_nz (or d_nnz) and o_nz (or o_nnz). By setting these parameters accurately, 3231 performance can be increased by more than a factor of 50. 3232 3233 Collective on Mat 3234 3235 Input Parameters: 3236 + B - the matrix 3237 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 3238 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 3239 . d_nz - number of block nonzeros per block row in diagonal portion of local 3240 submatrix (same for all local rows) 3241 . d_nnz - array containing the number of block nonzeros in the various block rows 3242 of the in diagonal portion of the local (possibly different for each block 3243 row) or NULL. If you plan to factor the matrix you must leave room for the diagonal entry and 3244 set it even if it is zero. 3245 . o_nz - number of block nonzeros per block row in the off-diagonal portion of local 3246 submatrix (same for all local rows). 3247 - o_nnz - array containing the number of nonzeros in the various block rows of the 3248 off-diagonal portion of the local submatrix (possibly different for 3249 each block row) or NULL. 3250 3251 If the *_nnz parameter is given then the *_nz parameter is ignored 3252 3253 Options Database Keys: 3254 + -mat_block_size - size of the blocks to use 3255 - -mat_use_hash_table <fact> 3256 3257 Notes: 3258 If PETSC_DECIDE or PETSC_DETERMINE is used for a particular argument on one processor 3259 than it must be used on all processors that share the object for that argument. 3260 3261 Storage Information: 3262 For a square global matrix we define each processor's diagonal portion 3263 to be its local rows and the corresponding columns (a square submatrix); 3264 each processor's off-diagonal portion encompasses the remainder of the 3265 local matrix (a rectangular submatrix). 3266 3267 The user can specify preallocated storage for the diagonal part of 3268 the local submatrix with either d_nz or d_nnz (not both). Set 3269 d_nz=PETSC_DEFAULT and d_nnz=NULL for PETSc to control dynamic 3270 memory allocation. Likewise, specify preallocated storage for the 3271 off-diagonal part of the local submatrix with o_nz or o_nnz (not both). 3272 3273 Consider a processor that owns rows 3, 4 and 5 of a parallel matrix. In 3274 the figure below we depict these three local rows and all columns (0-11). 3275 3276 .vb 3277 0 1 2 3 4 5 6 7 8 9 10 11 3278 -------------------------- 3279 row 3 |o o o d d d o o o o o o 3280 row 4 |o o o d d d o o o o o o 3281 row 5 |o o o d d d o o o o o o 3282 -------------------------- 3283 .ve 3284 3285 Thus, any entries in the d locations are stored in the d (diagonal) 3286 submatrix, and any entries in the o locations are stored in the 3287 o (off-diagonal) submatrix. Note that the d and the o submatrices are 3288 stored simply in the MATSEQBAIJ format for compressed row storage. 3289 3290 Now d_nz should indicate the number of block nonzeros per row in the d matrix, 3291 and o_nz should indicate the number of block nonzeros per row in the o matrix. 3292 In general, for PDE problems in which most nonzeros are near the diagonal, 3293 one expects d_nz >> o_nz. For large problems you MUST preallocate memory 3294 or you will get TERRIBLE performance; see the users' manual chapter on 3295 matrices. 3296 3297 You can call MatGetInfo() to get information on how effective the preallocation was; 3298 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3299 You can also run with the option -info and look for messages with the string 3300 malloc in them to see if additional memory allocation was needed. 3301 3302 Level: intermediate 3303 3304 .keywords: matrix, block, aij, compressed row, sparse, parallel 3305 3306 .seealso: MatCreate(), MatCreateSeqBAIJ(), MatSetValues(), MatCreateBAIJ(), MatMPIBAIJSetPreallocationCSR(), PetscSplitOwnership() 3307 @*/ 3308 PetscErrorCode MatMPIBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 3309 { 3310 PetscErrorCode ierr; 3311 3312 PetscFunctionBegin; 3313 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3314 PetscValidType(B,1); 3315 PetscValidLogicalCollectiveInt(B,bs,2); 3316 ierr = PetscTryMethod(B,"MatMPIBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,bs,d_nz,d_nnz,o_nz,o_nnz));CHKERRQ(ierr); 3317 PetscFunctionReturn(0); 3318 } 3319 3320 #undef __FUNCT__ 3321 #define __FUNCT__ "MatCreateBAIJ" 3322 /*@C 3323 MatCreateBAIJ - Creates a sparse parallel matrix in block AIJ format 3324 (block compressed row). For good matrix assembly performance 3325 the user should preallocate the matrix storage by setting the parameters 3326 d_nz (or d_nnz) and o_nz (or o_nnz). By setting these parameters accurately, 3327 performance can be increased by more than a factor of 50. 3328 3329 Collective on MPI_Comm 3330 3331 Input Parameters: 3332 + comm - MPI communicator 3333 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 3334 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 3335 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 3336 This value should be the same as the local size used in creating the 3337 y vector for the matrix-vector product y = Ax. 3338 . n - number of local columns (or PETSC_DECIDE to have calculated if N is given) 3339 This value should be the same as the local size used in creating the 3340 x vector for the matrix-vector product y = Ax. 3341 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 3342 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 3343 . d_nz - number of nonzero blocks per block row in diagonal portion of local 3344 submatrix (same for all local rows) 3345 . d_nnz - array containing the number of nonzero blocks in the various block rows 3346 of the in diagonal portion of the local (possibly different for each block 3347 row) or NULL. If you plan to factor the matrix you must leave room for the diagonal entry 3348 and set it even if it is zero. 3349 . o_nz - number of nonzero blocks per block row in the off-diagonal portion of local 3350 submatrix (same for all local rows). 3351 - o_nnz - array containing the number of nonzero blocks in the various block rows of the 3352 off-diagonal portion of the local submatrix (possibly different for 3353 each block row) or NULL. 3354 3355 Output Parameter: 3356 . A - the matrix 3357 3358 Options Database Keys: 3359 + -mat_block_size - size of the blocks to use 3360 - -mat_use_hash_table <fact> 3361 3362 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3363 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3364 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3365 3366 Notes: 3367 If the *_nnz parameter is given then the *_nz parameter is ignored 3368 3369 A nonzero block is any block that as 1 or more nonzeros in it 3370 3371 The user MUST specify either the local or global matrix dimensions 3372 (possibly both). 3373 3374 If PETSC_DECIDE or PETSC_DETERMINE is used for a particular argument on one processor 3375 than it must be used on all processors that share the object for that argument. 3376 3377 Storage Information: 3378 For a square global matrix we define each processor's diagonal portion 3379 to be its local rows and the corresponding columns (a square submatrix); 3380 each processor's off-diagonal portion encompasses the remainder of the 3381 local matrix (a rectangular submatrix). 3382 3383 The user can specify preallocated storage for the diagonal part of 3384 the local submatrix with either d_nz or d_nnz (not both). Set 3385 d_nz=PETSC_DEFAULT and d_nnz=NULL for PETSc to control dynamic 3386 memory allocation. Likewise, specify preallocated storage for the 3387 off-diagonal part of the local submatrix with o_nz or o_nnz (not both). 3388 3389 Consider a processor that owns rows 3, 4 and 5 of a parallel matrix. In 3390 the figure below we depict these three local rows and all columns (0-11). 3391 3392 .vb 3393 0 1 2 3 4 5 6 7 8 9 10 11 3394 -------------------------- 3395 row 3 |o o o d d d o o o o o o 3396 row 4 |o o o d d d o o o o o o 3397 row 5 |o o o d d d o o o o o o 3398 -------------------------- 3399 .ve 3400 3401 Thus, any entries in the d locations are stored in the d (diagonal) 3402 submatrix, and any entries in the o locations are stored in the 3403 o (off-diagonal) submatrix. Note that the d and the o submatrices are 3404 stored simply in the MATSEQBAIJ format for compressed row storage. 3405 3406 Now d_nz should indicate the number of block nonzeros per row in the d matrix, 3407 and o_nz should indicate the number of block nonzeros per row in the o matrix. 3408 In general, for PDE problems in which most nonzeros are near the diagonal, 3409 one expects d_nz >> o_nz. For large problems you MUST preallocate memory 3410 or you will get TERRIBLE performance; see the users' manual chapter on 3411 matrices. 3412 3413 Level: intermediate 3414 3415 .keywords: matrix, block, aij, compressed row, sparse, parallel 3416 3417 .seealso: MatCreate(), MatCreateSeqBAIJ(), MatSetValues(), MatCreateBAIJ(), MatMPIBAIJSetPreallocation(), MatMPIBAIJSetPreallocationCSR() 3418 @*/ 3419 PetscErrorCode MatCreateBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[],Mat *A) 3420 { 3421 PetscErrorCode ierr; 3422 PetscMPIInt size; 3423 3424 PetscFunctionBegin; 3425 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3426 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 3427 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3428 if (size > 1) { 3429 ierr = MatSetType(*A,MATMPIBAIJ);CHKERRQ(ierr); 3430 ierr = MatMPIBAIJSetPreallocation(*A,bs,d_nz,d_nnz,o_nz,o_nnz);CHKERRQ(ierr); 3431 } else { 3432 ierr = MatSetType(*A,MATSEQBAIJ);CHKERRQ(ierr); 3433 ierr = MatSeqBAIJSetPreallocation(*A,bs,d_nz,d_nnz);CHKERRQ(ierr); 3434 } 3435 PetscFunctionReturn(0); 3436 } 3437 3438 #undef __FUNCT__ 3439 #define __FUNCT__ "MatDuplicate_MPIBAIJ" 3440 static PetscErrorCode MatDuplicate_MPIBAIJ(Mat matin,MatDuplicateOption cpvalues,Mat *newmat) 3441 { 3442 Mat mat; 3443 Mat_MPIBAIJ *a,*oldmat = (Mat_MPIBAIJ*)matin->data; 3444 PetscErrorCode ierr; 3445 PetscInt len=0; 3446 3447 PetscFunctionBegin; 3448 *newmat = 0; 3449 ierr = MatCreate(PetscObjectComm((PetscObject)matin),&mat);CHKERRQ(ierr); 3450 ierr = MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N);CHKERRQ(ierr); 3451 ierr = MatSetType(mat,((PetscObject)matin)->type_name);CHKERRQ(ierr); 3452 ierr = PetscMemcpy(mat->ops,matin->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 3453 3454 mat->factortype = matin->factortype; 3455 mat->preallocated = PETSC_TRUE; 3456 mat->assembled = PETSC_TRUE; 3457 mat->insertmode = NOT_SET_VALUES; 3458 3459 a = (Mat_MPIBAIJ*)mat->data; 3460 mat->rmap->bs = matin->rmap->bs; 3461 a->bs2 = oldmat->bs2; 3462 a->mbs = oldmat->mbs; 3463 a->nbs = oldmat->nbs; 3464 a->Mbs = oldmat->Mbs; 3465 a->Nbs = oldmat->Nbs; 3466 3467 ierr = PetscLayoutReference(matin->rmap,&mat->rmap);CHKERRQ(ierr); 3468 ierr = PetscLayoutReference(matin->cmap,&mat->cmap);CHKERRQ(ierr); 3469 3470 a->size = oldmat->size; 3471 a->rank = oldmat->rank; 3472 a->donotstash = oldmat->donotstash; 3473 a->roworiented = oldmat->roworiented; 3474 a->rowindices = 0; 3475 a->rowvalues = 0; 3476 a->getrowactive = PETSC_FALSE; 3477 a->barray = 0; 3478 a->rstartbs = oldmat->rstartbs; 3479 a->rendbs = oldmat->rendbs; 3480 a->cstartbs = oldmat->cstartbs; 3481 a->cendbs = oldmat->cendbs; 3482 3483 /* hash table stuff */ 3484 a->ht = 0; 3485 a->hd = 0; 3486 a->ht_size = 0; 3487 a->ht_flag = oldmat->ht_flag; 3488 a->ht_fact = oldmat->ht_fact; 3489 a->ht_total_ct = 0; 3490 a->ht_insert_ct = 0; 3491 3492 ierr = PetscMemcpy(a->rangebs,oldmat->rangebs,(a->size+1)*sizeof(PetscInt));CHKERRQ(ierr); 3493 if (oldmat->colmap) { 3494 #if defined(PETSC_USE_CTABLE) 3495 ierr = PetscTableCreateCopy(oldmat->colmap,&a->colmap);CHKERRQ(ierr); 3496 #else 3497 ierr = PetscMalloc1(a->Nbs,&a->colmap);CHKERRQ(ierr); 3498 ierr = PetscLogObjectMemory((PetscObject)mat,(a->Nbs)*sizeof(PetscInt));CHKERRQ(ierr); 3499 ierr = PetscMemcpy(a->colmap,oldmat->colmap,(a->Nbs)*sizeof(PetscInt));CHKERRQ(ierr); 3500 #endif 3501 } else a->colmap = 0; 3502 3503 if (oldmat->garray && (len = ((Mat_SeqBAIJ*)(oldmat->B->data))->nbs)) { 3504 ierr = PetscMalloc1(len,&a->garray);CHKERRQ(ierr); 3505 ierr = PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt));CHKERRQ(ierr); 3506 ierr = PetscMemcpy(a->garray,oldmat->garray,len*sizeof(PetscInt));CHKERRQ(ierr); 3507 } else a->garray = 0; 3508 3509 ierr = MatStashCreate_Private(PetscObjectComm((PetscObject)matin),matin->rmap->bs,&mat->bstash);CHKERRQ(ierr); 3510 ierr = VecDuplicate(oldmat->lvec,&a->lvec);CHKERRQ(ierr); 3511 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec);CHKERRQ(ierr); 3512 ierr = VecScatterCopy(oldmat->Mvctx,&a->Mvctx);CHKERRQ(ierr); 3513 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx);CHKERRQ(ierr); 3514 3515 ierr = MatDuplicate(oldmat->A,cpvalues,&a->A);CHKERRQ(ierr); 3516 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A);CHKERRQ(ierr); 3517 ierr = MatDuplicate(oldmat->B,cpvalues,&a->B);CHKERRQ(ierr); 3518 ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B);CHKERRQ(ierr); 3519 ierr = PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist);CHKERRQ(ierr); 3520 *newmat = mat; 3521 PetscFunctionReturn(0); 3522 } 3523 3524 #undef __FUNCT__ 3525 #define __FUNCT__ "MatLoad_MPIBAIJ" 3526 PetscErrorCode MatLoad_MPIBAIJ(Mat newmat,PetscViewer viewer) 3527 { 3528 PetscErrorCode ierr; 3529 int fd; 3530 PetscInt i,nz,j,rstart,rend; 3531 PetscScalar *vals,*buf; 3532 MPI_Comm comm; 3533 MPI_Status status; 3534 PetscMPIInt rank,size,maxnz; 3535 PetscInt header[4],*rowlengths = 0,M,N,m,*rowners,*cols; 3536 PetscInt *locrowlens = NULL,*procsnz = NULL,*browners = NULL; 3537 PetscInt jj,*mycols,*ibuf,bs = newmat->rmap->bs,Mbs,mbs,extra_rows,mmax; 3538 PetscMPIInt tag = ((PetscObject)viewer)->tag; 3539 PetscInt *dlens = NULL,*odlens = NULL,*mask = NULL,*masked1 = NULL,*masked2 = NULL,rowcount,odcount; 3540 PetscInt dcount,kmax,k,nzcount,tmp,mend; 3541 3542 PetscFunctionBegin; 3543 /* force binary viewer to load .info file if it has not yet done so */ 3544 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 3545 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3546 ierr = PetscOptionsBegin(comm,NULL,"Options for loading MPIBAIJ matrix 2","Mat");CHKERRQ(ierr); 3547 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 3548 ierr = PetscOptionsEnd();CHKERRQ(ierr); 3549 if (bs < 0) bs = 1; 3550 3551 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3552 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 3553 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3554 if (!rank) { 3555 ierr = PetscBinaryRead(fd,(char*)header,4,PETSC_INT);CHKERRQ(ierr); 3556 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object"); 3557 if (header[3] < 0) SETERRQ(PetscObjectComm((PetscObject)newmat),PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk, cannot load as MPIAIJ"); 3558 } 3559 ierr = MPI_Bcast(header+1,3,MPIU_INT,0,comm);CHKERRQ(ierr); 3560 M = header[1]; N = header[2]; 3561 3562 /* If global sizes are set, check if they are consistent with that given in the file */ 3563 if (newmat->rmap->N >= 0 && newmat->rmap->N != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Inconsistent # of rows:Matrix in file has (%D) and input matrix has (%D)",newmat->rmap->N,M); 3564 if (newmat->cmap->N >= 0 && newmat->cmap->N != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Inconsistent # of cols:Matrix in file has (%D) and input matrix has (%D)",newmat->cmap->N,N); 3565 3566 if (M != N) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Can only do square matrices"); 3567 3568 /* 3569 This code adds extra rows to make sure the number of rows is 3570 divisible by the blocksize 3571 */ 3572 Mbs = M/bs; 3573 extra_rows = bs - M + bs*Mbs; 3574 if (extra_rows == bs) extra_rows = 0; 3575 else Mbs++; 3576 if (extra_rows && !rank) { 3577 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 3578 } 3579 3580 /* determine ownership of all rows */ 3581 if (newmat->rmap->n < 0) { /* PETSC_DECIDE */ 3582 mbs = Mbs/size + ((Mbs % size) > rank); 3583 m = mbs*bs; 3584 } else { /* User set */ 3585 m = newmat->rmap->n; 3586 mbs = m/bs; 3587 } 3588 ierr = PetscMalloc2(size+1,&rowners,size+1,&browners);CHKERRQ(ierr); 3589 ierr = MPI_Allgather(&mbs,1,MPIU_INT,rowners+1,1,MPIU_INT,comm);CHKERRQ(ierr); 3590 3591 /* process 0 needs enough room for process with most rows */ 3592 if (!rank) { 3593 mmax = rowners[1]; 3594 for (i=2; i<=size; i++) { 3595 mmax = PetscMax(mmax,rowners[i]); 3596 } 3597 mmax*=bs; 3598 } else mmax = -1; /* unused, but compiler warns anyway */ 3599 3600 rowners[0] = 0; 3601 for (i=2; i<=size; i++) rowners[i] += rowners[i-1]; 3602 for (i=0; i<=size; i++) browners[i] = rowners[i]*bs; 3603 rstart = rowners[rank]; 3604 rend = rowners[rank+1]; 3605 3606 /* distribute row lengths to all processors */ 3607 ierr = PetscMalloc1(m,&locrowlens);CHKERRQ(ierr); 3608 if (!rank) { 3609 mend = m; 3610 if (size == 1) mend = mend - extra_rows; 3611 ierr = PetscBinaryRead(fd,locrowlens,mend,PETSC_INT);CHKERRQ(ierr); 3612 for (j=mend; j<m; j++) locrowlens[j] = 1; 3613 ierr = PetscMalloc1(mmax,&rowlengths);CHKERRQ(ierr); 3614 ierr = PetscCalloc1(size,&procsnz);CHKERRQ(ierr); 3615 for (j=0; j<m; j++) { 3616 procsnz[0] += locrowlens[j]; 3617 } 3618 for (i=1; i<size; i++) { 3619 mend = browners[i+1] - browners[i]; 3620 if (i == size-1) mend = mend - extra_rows; 3621 ierr = PetscBinaryRead(fd,rowlengths,mend,PETSC_INT);CHKERRQ(ierr); 3622 for (j=mend; j<browners[i+1] - browners[i]; j++) rowlengths[j] = 1; 3623 /* calculate the number of nonzeros on each processor */ 3624 for (j=0; j<browners[i+1]-browners[i]; j++) { 3625 procsnz[i] += rowlengths[j]; 3626 } 3627 ierr = MPI_Send(rowlengths,browners[i+1]-browners[i],MPIU_INT,i,tag,comm);CHKERRQ(ierr); 3628 } 3629 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3630 } else { 3631 ierr = MPI_Recv(locrowlens,m,MPIU_INT,0,tag,comm,&status);CHKERRQ(ierr); 3632 } 3633 3634 if (!rank) { 3635 /* determine max buffer needed and allocate it */ 3636 maxnz = procsnz[0]; 3637 for (i=1; i<size; i++) { 3638 maxnz = PetscMax(maxnz,procsnz[i]); 3639 } 3640 ierr = PetscMalloc1(maxnz,&cols);CHKERRQ(ierr); 3641 3642 /* read in my part of the matrix column indices */ 3643 nz = procsnz[0]; 3644 ierr = PetscMalloc1(nz+1,&ibuf);CHKERRQ(ierr); 3645 mycols = ibuf; 3646 if (size == 1) nz -= extra_rows; 3647 ierr = PetscBinaryRead(fd,mycols,nz,PETSC_INT);CHKERRQ(ierr); 3648 if (size == 1) { 3649 for (i=0; i< extra_rows; i++) mycols[nz+i] = M+i; 3650 } 3651 3652 /* read in every ones (except the last) and ship off */ 3653 for (i=1; i<size-1; i++) { 3654 nz = procsnz[i]; 3655 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 3656 ierr = MPI_Send(cols,nz,MPIU_INT,i,tag,comm);CHKERRQ(ierr); 3657 } 3658 /* read in the stuff for the last proc */ 3659 if (size != 1) { 3660 nz = procsnz[size-1] - extra_rows; /* the extra rows are not on the disk */ 3661 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 3662 for (i=0; i<extra_rows; i++) cols[nz+i] = M+i; 3663 ierr = MPI_Send(cols,nz+extra_rows,MPIU_INT,size-1,tag,comm);CHKERRQ(ierr); 3664 } 3665 ierr = PetscFree(cols);CHKERRQ(ierr); 3666 } else { 3667 /* determine buffer space needed for message */ 3668 nz = 0; 3669 for (i=0; i<m; i++) { 3670 nz += locrowlens[i]; 3671 } 3672 ierr = PetscMalloc1(nz+1,&ibuf);CHKERRQ(ierr); 3673 mycols = ibuf; 3674 /* receive message of column indices*/ 3675 ierr = MPI_Recv(mycols,nz,MPIU_INT,0,tag,comm,&status);CHKERRQ(ierr); 3676 ierr = MPI_Get_count(&status,MPIU_INT,&maxnz);CHKERRQ(ierr); 3677 if (maxnz != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"something is wrong with file"); 3678 } 3679 3680 /* loop over local rows, determining number of off diagonal entries */ 3681 ierr = PetscMalloc2(rend-rstart,&dlens,rend-rstart,&odlens);CHKERRQ(ierr); 3682 ierr = PetscCalloc3(Mbs,&mask,Mbs,&masked1,Mbs,&masked2);CHKERRQ(ierr); 3683 rowcount = 0; nzcount = 0; 3684 for (i=0; i<mbs; i++) { 3685 dcount = 0; 3686 odcount = 0; 3687 for (j=0; j<bs; j++) { 3688 kmax = locrowlens[rowcount]; 3689 for (k=0; k<kmax; k++) { 3690 tmp = mycols[nzcount++]/bs; 3691 if (!mask[tmp]) { 3692 mask[tmp] = 1; 3693 if (tmp < rstart || tmp >= rend) masked2[odcount++] = tmp; 3694 else masked1[dcount++] = tmp; 3695 } 3696 } 3697 rowcount++; 3698 } 3699 3700 dlens[i] = dcount; 3701 odlens[i] = odcount; 3702 3703 /* zero out the mask elements we set */ 3704 for (j=0; j<dcount; j++) mask[masked1[j]] = 0; 3705 for (j=0; j<odcount; j++) mask[masked2[j]] = 0; 3706 } 3707 3708 ierr = MatSetSizes(newmat,m,m,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 3709 ierr = MatMPIBAIJSetPreallocation(newmat,bs,0,dlens,0,odlens);CHKERRQ(ierr); 3710 3711 if (!rank) { 3712 ierr = PetscMalloc1(maxnz+1,&buf);CHKERRQ(ierr); 3713 /* read in my part of the matrix numerical values */ 3714 nz = procsnz[0]; 3715 vals = buf; 3716 mycols = ibuf; 3717 if (size == 1) nz -= extra_rows; 3718 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 3719 if (size == 1) { 3720 for (i=0; i< extra_rows; i++) vals[nz+i] = 1.0; 3721 } 3722 3723 /* insert into matrix */ 3724 jj = rstart*bs; 3725 for (i=0; i<m; i++) { 3726 ierr = MatSetValues_MPIBAIJ(newmat,1,&jj,locrowlens[i],mycols,vals,INSERT_VALUES);CHKERRQ(ierr); 3727 mycols += locrowlens[i]; 3728 vals += locrowlens[i]; 3729 jj++; 3730 } 3731 /* read in other processors (except the last one) and ship out */ 3732 for (i=1; i<size-1; i++) { 3733 nz = procsnz[i]; 3734 vals = buf; 3735 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 3736 ierr = MPIULong_Send(vals,nz,MPIU_SCALAR,i,((PetscObject)newmat)->tag,comm);CHKERRQ(ierr); 3737 } 3738 /* the last proc */ 3739 if (size != 1) { 3740 nz = procsnz[i] - extra_rows; 3741 vals = buf; 3742 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 3743 for (i=0; i<extra_rows; i++) vals[nz+i] = 1.0; 3744 ierr = MPIULong_Send(vals,nz+extra_rows,MPIU_SCALAR,size-1,((PetscObject)newmat)->tag,comm);CHKERRQ(ierr); 3745 } 3746 ierr = PetscFree(procsnz);CHKERRQ(ierr); 3747 } else { 3748 /* receive numeric values */ 3749 ierr = PetscMalloc1(nz+1,&buf);CHKERRQ(ierr); 3750 3751 /* receive message of values*/ 3752 vals = buf; 3753 mycols = ibuf; 3754 ierr = MPIULong_Recv(vals,nz,MPIU_SCALAR,0,((PetscObject)newmat)->tag,comm);CHKERRQ(ierr); 3755 3756 /* insert into matrix */ 3757 jj = rstart*bs; 3758 for (i=0; i<m; i++) { 3759 ierr = MatSetValues_MPIBAIJ(newmat,1,&jj,locrowlens[i],mycols,vals,INSERT_VALUES);CHKERRQ(ierr); 3760 mycols += locrowlens[i]; 3761 vals += locrowlens[i]; 3762 jj++; 3763 } 3764 } 3765 ierr = PetscFree(locrowlens);CHKERRQ(ierr); 3766 ierr = PetscFree(buf);CHKERRQ(ierr); 3767 ierr = PetscFree(ibuf);CHKERRQ(ierr); 3768 ierr = PetscFree2(rowners,browners);CHKERRQ(ierr); 3769 ierr = PetscFree2(dlens,odlens);CHKERRQ(ierr); 3770 ierr = PetscFree3(mask,masked1,masked2);CHKERRQ(ierr); 3771 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3772 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3773 PetscFunctionReturn(0); 3774 } 3775 3776 #undef __FUNCT__ 3777 #define __FUNCT__ "MatMPIBAIJSetHashTableFactor" 3778 /*@ 3779 MatMPIBAIJSetHashTableFactor - Sets the factor required to compute the size of the HashTable. 3780 3781 Input Parameters: 3782 . mat - the matrix 3783 . fact - factor 3784 3785 Not Collective, each process can use a different factor 3786 3787 Level: advanced 3788 3789 Notes: 3790 This can also be set by the command line option: -mat_use_hash_table <fact> 3791 3792 .keywords: matrix, hashtable, factor, HT 3793 3794 .seealso: MatSetOption() 3795 @*/ 3796 PetscErrorCode MatMPIBAIJSetHashTableFactor(Mat mat,PetscReal fact) 3797 { 3798 PetscErrorCode ierr; 3799 3800 PetscFunctionBegin; 3801 ierr = PetscTryMethod(mat,"MatSetHashTableFactor_C",(Mat,PetscReal),(mat,fact));CHKERRQ(ierr); 3802 PetscFunctionReturn(0); 3803 } 3804 3805 #undef __FUNCT__ 3806 #define __FUNCT__ "MatSetHashTableFactor_MPIBAIJ" 3807 PetscErrorCode MatSetHashTableFactor_MPIBAIJ(Mat mat,PetscReal fact) 3808 { 3809 Mat_MPIBAIJ *baij; 3810 3811 PetscFunctionBegin; 3812 baij = (Mat_MPIBAIJ*)mat->data; 3813 baij->ht_fact = fact; 3814 PetscFunctionReturn(0); 3815 } 3816 3817 #undef __FUNCT__ 3818 #define __FUNCT__ "MatMPIBAIJGetSeqBAIJ" 3819 PetscErrorCode MatMPIBAIJGetSeqBAIJ(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[]) 3820 { 3821 Mat_MPIBAIJ *a = (Mat_MPIBAIJ*)A->data; 3822 3823 PetscFunctionBegin; 3824 if (Ad) *Ad = a->A; 3825 if (Ao) *Ao = a->B; 3826 if (colmap) *colmap = a->garray; 3827 PetscFunctionReturn(0); 3828 } 3829 3830 /* 3831 Special version for direct calls from Fortran (to eliminate two function call overheads 3832 */ 3833 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3834 #define matmpibaijsetvaluesblocked_ MATMPIBAIJSETVALUESBLOCKED 3835 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3836 #define matmpibaijsetvaluesblocked_ matmpibaijsetvaluesblocked 3837 #endif 3838 3839 #undef __FUNCT__ 3840 #define __FUNCT__ "matmpibiajsetvaluesblocked" 3841 /*@C 3842 MatMPIBAIJSetValuesBlocked - Direct Fortran call to replace call to MatSetValuesBlocked() 3843 3844 Collective on Mat 3845 3846 Input Parameters: 3847 + mat - the matrix 3848 . min - number of input rows 3849 . im - input rows 3850 . nin - number of input columns 3851 . in - input columns 3852 . v - numerical values input 3853 - addvin - INSERT_VALUES or ADD_VALUES 3854 3855 Notes: This has a complete copy of MatSetValuesBlocked_MPIBAIJ() which is terrible code un-reuse. 3856 3857 Level: advanced 3858 3859 .seealso: MatSetValuesBlocked() 3860 @*/ 3861 PetscErrorCode matmpibaijsetvaluesblocked_(Mat *matin,PetscInt *min,const PetscInt im[],PetscInt *nin,const PetscInt in[],const MatScalar v[],InsertMode *addvin) 3862 { 3863 /* convert input arguments to C version */ 3864 Mat mat = *matin; 3865 PetscInt m = *min, n = *nin; 3866 InsertMode addv = *addvin; 3867 3868 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat->data; 3869 const MatScalar *value; 3870 MatScalar *barray = baij->barray; 3871 PetscBool roworiented = baij->roworiented; 3872 PetscErrorCode ierr; 3873 PetscInt i,j,ii,jj,row,col,rstart=baij->rstartbs; 3874 PetscInt rend=baij->rendbs,cstart=baij->cstartbs,stepval; 3875 PetscInt cend=baij->cendbs,bs=mat->rmap->bs,bs2=baij->bs2; 3876 3877 PetscFunctionBegin; 3878 /* tasks normally handled by MatSetValuesBlocked() */ 3879 if (mat->insertmode == NOT_SET_VALUES) mat->insertmode = addv; 3880 #if defined(PETSC_USE_DEBUG) 3881 else if (mat->insertmode != addv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values"); 3882 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3883 #endif 3884 if (mat->assembled) { 3885 mat->was_assembled = PETSC_TRUE; 3886 mat->assembled = PETSC_FALSE; 3887 } 3888 ierr = PetscLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 3889 3890 3891 if (!barray) { 3892 ierr = PetscMalloc1(bs2,&barray);CHKERRQ(ierr); 3893 baij->barray = barray; 3894 } 3895 3896 if (roworiented) stepval = (n-1)*bs; 3897 else stepval = (m-1)*bs; 3898 3899 for (i=0; i<m; i++) { 3900 if (im[i] < 0) continue; 3901 #if defined(PETSC_USE_DEBUG) 3902 if (im[i] >= baij->Mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large, row %D max %D",im[i],baij->Mbs-1); 3903 #endif 3904 if (im[i] >= rstart && im[i] < rend) { 3905 row = im[i] - rstart; 3906 for (j=0; j<n; j++) { 3907 /* If NumCol = 1 then a copy is not required */ 3908 if ((roworiented) && (n == 1)) { 3909 barray = (MatScalar*)v + i*bs2; 3910 } else if ((!roworiented) && (m == 1)) { 3911 barray = (MatScalar*)v + j*bs2; 3912 } else { /* Here a copy is required */ 3913 if (roworiented) { 3914 value = v + i*(stepval+bs)*bs + j*bs; 3915 } else { 3916 value = v + j*(stepval+bs)*bs + i*bs; 3917 } 3918 for (ii=0; ii<bs; ii++,value+=stepval) { 3919 for (jj=0; jj<bs; jj++) { 3920 *barray++ = *value++; 3921 } 3922 } 3923 barray -=bs2; 3924 } 3925 3926 if (in[j] >= cstart && in[j] < cend) { 3927 col = in[j] - cstart; 3928 ierr = MatSetValuesBlocked_SeqBAIJ_Inlined(baij->A,row,col,barray,addv,im[i],in[j]);CHKERRQ(ierr); 3929 } else if (in[j] < 0) continue; 3930 #if defined(PETSC_USE_DEBUG) 3931 else if (in[j] >= baij->Nbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large, col %D max %D",in[j],baij->Nbs-1); 3932 #endif 3933 else { 3934 if (mat->was_assembled) { 3935 if (!baij->colmap) { 3936 ierr = MatCreateColmap_MPIBAIJ_Private(mat);CHKERRQ(ierr); 3937 } 3938 3939 #if defined(PETSC_USE_DEBUG) 3940 #if defined(PETSC_USE_CTABLE) 3941 { PetscInt data; 3942 ierr = PetscTableFind(baij->colmap,in[j]+1,&data);CHKERRQ(ierr); 3943 if ((data - 1) % bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Incorrect colmap"); 3944 } 3945 #else 3946 if ((baij->colmap[in[j]] - 1) % bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Incorrect colmap"); 3947 #endif 3948 #endif 3949 #if defined(PETSC_USE_CTABLE) 3950 ierr = PetscTableFind(baij->colmap,in[j]+1,&col);CHKERRQ(ierr); 3951 col = (col - 1)/bs; 3952 #else 3953 col = (baij->colmap[in[j]] - 1)/bs; 3954 #endif 3955 if (col < 0 && !((Mat_SeqBAIJ*)(baij->A->data))->nonew) { 3956 ierr = MatDisAssemble_MPIBAIJ(mat);CHKERRQ(ierr); 3957 col = in[j]; 3958 } 3959 } else col = in[j]; 3960 ierr = MatSetValuesBlocked_SeqBAIJ_Inlined(baij->B,row,col,barray,addv,im[i],in[j]);CHKERRQ(ierr); 3961 } 3962 } 3963 } else { 3964 if (!baij->donotstash) { 3965 if (roworiented) { 3966 ierr = MatStashValuesRowBlocked_Private(&mat->bstash,im[i],n,in,v,m,n,i);CHKERRQ(ierr); 3967 } else { 3968 ierr = MatStashValuesColBlocked_Private(&mat->bstash,im[i],n,in,v,m,n,i);CHKERRQ(ierr); 3969 } 3970 } 3971 } 3972 } 3973 3974 /* task normally handled by MatSetValuesBlocked() */ 3975 ierr = PetscLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 3976 PetscFunctionReturn(0); 3977 } 3978 3979 #undef __FUNCT__ 3980 #define __FUNCT__ "MatCreateMPIBAIJWithArrays" 3981 /*@ 3982 MatCreateMPIBAIJWithArrays - creates a MPI BAIJ matrix using arrays that contain in standard 3983 CSR format the local rows. 3984 3985 Collective on MPI_Comm 3986 3987 Input Parameters: 3988 + comm - MPI communicator 3989 . bs - the block size, only a block size of 1 is supported 3990 . m - number of local rows (Cannot be PETSC_DECIDE) 3991 . n - This value should be the same as the local size used in creating the 3992 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 3993 calculated if N is given) For square matrices n is almost always m. 3994 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 3995 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 3996 . i - row indices 3997 . j - column indices 3998 - a - matrix values 3999 4000 Output Parameter: 4001 . mat - the matrix 4002 4003 Level: intermediate 4004 4005 Notes: 4006 The i, j, and a arrays ARE copied by this routine into the internal format used by PETSc; 4007 thus you CANNOT change the matrix entries by changing the values of a[] after you have 4008 called this routine. Use MatCreateMPIAIJWithSplitArrays() to avoid needing to copy the arrays. 4009 4010 The order of the entries in values is the same as the block compressed sparse row storage format; that is, it is 4011 the same as a three dimensional array in Fortran values(bs,bs,nnz) that contains the first column of the first 4012 block, followed by the second column of the first block etc etc. That is, the blocks are contiguous in memory 4013 with column-major ordering within blocks. 4014 4015 The i and j indices are 0 based, and i indices are indices corresponding to the local j array. 4016 4017 .keywords: matrix, aij, compressed row, sparse, parallel 4018 4019 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatMPIAIJSetPreallocation(), MatMPIAIJSetPreallocationCSR(), 4020 MPIAIJ, MatCreateAIJ(), MatCreateMPIAIJWithSplitArrays() 4021 @*/ 4022 PetscErrorCode MatCreateMPIBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,const PetscInt i[],const PetscInt j[],const PetscScalar a[],Mat *mat) 4023 { 4024 PetscErrorCode ierr; 4025 4026 PetscFunctionBegin; 4027 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4028 if (m < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"local number of rows (m) cannot be PETSC_DECIDE, or negative"); 4029 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4030 ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr); 4031 ierr = MatSetType(*mat,MATMPISBAIJ);CHKERRQ(ierr); 4032 ierr = MatSetOption(*mat,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 4033 ierr = MatMPIBAIJSetPreallocationCSR(*mat,bs,i,j,a);CHKERRQ(ierr); 4034 ierr = MatSetOption(*mat,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 4035 PetscFunctionReturn(0); 4036 } 4037 4038 #undef __FUNCT__ 4039 #define __FUNCT__ "MatCreateMPIMatConcatenateSeqMat_MPIBAIJ" 4040 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_MPIBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4041 { 4042 PetscErrorCode ierr; 4043 PetscInt m,N,i,rstart,nnz,Ii,bs,cbs; 4044 PetscInt *indx; 4045 PetscScalar *values; 4046 4047 PetscFunctionBegin; 4048 ierr = MatGetSize(inmat,&m,&N);CHKERRQ(ierr); 4049 if (scall == MAT_INITIAL_MATRIX) { /* symbolic phase */ 4050 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)inmat->data; 4051 PetscInt *dnz,*onz,sum,mbs,Nbs; 4052 PetscInt *bindx,rmax=a->rmax,j; 4053 4054 ierr = MatGetBlockSizes(inmat,&bs,&cbs);CHKERRQ(ierr); 4055 mbs = m/bs; Nbs = N/cbs; 4056 if (n == PETSC_DECIDE) { 4057 ierr = PetscSplitOwnership(comm,&n,&Nbs);CHKERRQ(ierr); 4058 } 4059 /* Check sum(n) = Nbs */ 4060 ierr = MPIU_Allreduce(&n,&sum,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); 4061 if (sum != Nbs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Sum of local columns != global columns %d",Nbs); 4062 4063 ierr = MPI_Scan(&mbs, &rstart,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); 4064 rstart -= mbs; 4065 4066 ierr = PetscMalloc1(rmax,&bindx);CHKERRQ(ierr); 4067 ierr = MatPreallocateInitialize(comm,mbs,n,dnz,onz);CHKERRQ(ierr); 4068 for (i=0; i<mbs; i++) { 4069 ierr = MatGetRow_SeqBAIJ(inmat,i*bs,&nnz,&indx,NULL);CHKERRQ(ierr); /* non-blocked nnz and indx */ 4070 nnz = nnz/bs; 4071 for (j=0; j<nnz; j++) bindx[j] = indx[j*bs]/bs; 4072 ierr = MatPreallocateSet(i+rstart,nnz,bindx,dnz,onz);CHKERRQ(ierr); 4073 ierr = MatRestoreRow_SeqBAIJ(inmat,i*bs,&nnz,&indx,NULL);CHKERRQ(ierr); 4074 } 4075 ierr = PetscFree(bindx);CHKERRQ(ierr); 4076 4077 ierr = MatCreate(comm,outmat);CHKERRQ(ierr); 4078 ierr = MatSetSizes(*outmat,m,n*bs,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 4079 ierr = MatSetBlockSizes(*outmat,bs,cbs);CHKERRQ(ierr); 4080 ierr = MatSetType(*outmat,MATMPIBAIJ);CHKERRQ(ierr); 4081 ierr = MatMPIBAIJSetPreallocation(*outmat,bs,0,dnz,0,onz);CHKERRQ(ierr); 4082 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 4083 } 4084 4085 /* numeric phase */ 4086 ierr = MatGetBlockSizes(inmat,&bs,&cbs);CHKERRQ(ierr); 4087 ierr = MatGetOwnershipRange(*outmat,&rstart,NULL);CHKERRQ(ierr); 4088 4089 for (i=0; i<m; i++) { 4090 ierr = MatGetRow_SeqBAIJ(inmat,i,&nnz,&indx,&values);CHKERRQ(ierr); 4091 Ii = i + rstart; 4092 ierr = MatSetValues(*outmat,1,&Ii,nnz,indx,values,INSERT_VALUES);CHKERRQ(ierr); 4093 ierr = MatRestoreRow_SeqBAIJ(inmat,i,&nnz,&indx,&values);CHKERRQ(ierr); 4094 } 4095 ierr = MatAssemblyBegin(*outmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4096 ierr = MatAssemblyEnd(*outmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4097 PetscFunctionReturn(0); 4098 } 4099