1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: baij.c,v 1.148 1998/12/03 04:01:06 bsmith Exp balay $"; 3 #endif 4 5 /* 6 Defines the basic matrix operations for the BAIJ (compressed row) 7 matrix storage format. 8 */ 9 #include "sys.h" 10 #include "src/mat/impls/baij/seq/baij.h" 11 #include "src/vec/vecimpl.h" 12 #include "src/inline/spops.h" 13 14 #define CHUNKSIZE 10 15 16 #undef __FUNC__ 17 #define __FUNC__ "MatMarkDiag_SeqBAIJ" 18 int MatMarkDiag_SeqBAIJ(Mat A) 19 { 20 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 21 int i,j, *diag, m = a->mbs; 22 23 PetscFunctionBegin; 24 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 25 PLogObjectMemory(A,(m+1)*sizeof(int)); 26 for ( i=0; i<m; i++ ) { 27 diag[i] = a->i[i+1]; 28 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 29 if (a->j[j] == i) { 30 diag[i] = j; 31 break; 32 } 33 } 34 } 35 a->diag = diag; 36 PetscFunctionReturn(0); 37 } 38 39 40 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 41 42 #undef __FUNC__ 43 #define __FUNC__ "MatGetRowIJ_SeqBAIJ" 44 static int MatGetRowIJ_SeqBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 45 PetscTruth *done) 46 { 47 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 48 int ierr, n = a->mbs,i; 49 50 PetscFunctionBegin; 51 *nn = n; 52 if (!ia) PetscFunctionReturn(0); 53 if (symmetric) { 54 ierr = MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 55 } else if (oshift == 1) { 56 /* temporarily add 1 to i and j indices */ 57 int nz = a->i[n] + 1; 58 for ( i=0; i<nz; i++ ) a->j[i]++; 59 for ( i=0; i<n+1; i++ ) a->i[i]++; 60 *ia = a->i; *ja = a->j; 61 } else { 62 *ia = a->i; *ja = a->j; 63 } 64 65 PetscFunctionReturn(0); 66 } 67 68 #undef __FUNC__ 69 #define __FUNC__ "MatRestoreRowIJ_SeqBAIJ" 70 static int MatRestoreRowIJ_SeqBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 71 PetscTruth *done) 72 { 73 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 74 int i,n = a->mbs; 75 76 PetscFunctionBegin; 77 if (!ia) PetscFunctionReturn(0); 78 if (symmetric) { 79 PetscFree(*ia); 80 PetscFree(*ja); 81 } else if (oshift == 1) { 82 int nz = a->i[n]; 83 for ( i=0; i<nz; i++ ) a->j[i]--; 84 for ( i=0; i<n+1; i++ ) a->i[i]--; 85 } 86 PetscFunctionReturn(0); 87 } 88 89 #undef __FUNC__ 90 #define __FUNC__ "MatGetBlockSize_SeqBAIJ" 91 int MatGetBlockSize_SeqBAIJ(Mat mat, int *bs) 92 { 93 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *) mat->data; 94 95 PetscFunctionBegin; 96 *bs = baij->bs; 97 PetscFunctionReturn(0); 98 } 99 100 101 #undef __FUNC__ 102 #define __FUNC__ "MatDestroy_SeqBAIJ" 103 int MatDestroy_SeqBAIJ(Mat A) 104 { 105 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 106 int ierr; 107 108 if (--A->refct > 0) PetscFunctionReturn(0); 109 110 if (A->mapping) { 111 ierr = ISLocalToGlobalMappingDestroy(A->mapping); CHKERRQ(ierr); 112 } 113 if (A->bmapping) { 114 ierr = ISLocalToGlobalMappingDestroy(A->bmapping); CHKERRQ(ierr); 115 } 116 if (A->rmap) { 117 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 118 } 119 if (A->cmap) { 120 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 121 } 122 #if defined(USE_PETSC_LOG) 123 PLogObjectState((PetscObject) A,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 124 #endif 125 PetscFree(a->a); 126 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 127 if (a->diag) PetscFree(a->diag); 128 if (a->ilen) PetscFree(a->ilen); 129 if (a->imax) PetscFree(a->imax); 130 if (a->solve_work) PetscFree(a->solve_work); 131 if (a->mult_work) PetscFree(a->mult_work); 132 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 133 PetscFree(a); 134 PLogObjectDestroy(A); 135 PetscHeaderDestroy(A); 136 PetscFunctionReturn(0); 137 } 138 139 #undef __FUNC__ 140 #define __FUNC__ "MatSetOption_SeqBAIJ" 141 int MatSetOption_SeqBAIJ(Mat A,MatOption op) 142 { 143 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 144 145 PetscFunctionBegin; 146 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 147 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 148 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 149 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 150 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 151 else if (op == MAT_NEW_NONZERO_LOCATION_ERROR) a->nonew = -1; 152 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERROR) a->nonew = -2; 153 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 154 else if (op == MAT_ROWS_SORTED || 155 op == MAT_ROWS_UNSORTED || 156 op == MAT_SYMMETRIC || 157 op == MAT_STRUCTURALLY_SYMMETRIC || 158 op == MAT_YES_NEW_DIAGONALS || 159 op == MAT_IGNORE_OFF_PROC_ENTRIES || 160 op == MAT_USE_HASH_TABLE) { 161 PLogInfo(A,"MatSetOption_SeqBAIJ:Option ignored\n"); 162 } else if (op == MAT_NO_NEW_DIAGONALS) { 163 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 164 } else { 165 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 166 } 167 PetscFunctionReturn(0); 168 } 169 170 171 #undef __FUNC__ 172 #define __FUNC__ "MatGetSize_SeqBAIJ" 173 int MatGetSize_SeqBAIJ(Mat A,int *m,int *n) 174 { 175 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 176 177 PetscFunctionBegin; 178 if (m) *m = a->m; 179 if (n) *n = a->n; 180 PetscFunctionReturn(0); 181 } 182 183 #undef __FUNC__ 184 #define __FUNC__ "MatGetOwnershipRange_SeqBAIJ" 185 int MatGetOwnershipRange_SeqBAIJ(Mat A,int *m,int *n) 186 { 187 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 188 189 PetscFunctionBegin; 190 *m = 0; *n = a->m; 191 PetscFunctionReturn(0); 192 } 193 194 #undef __FUNC__ 195 #define __FUNC__ "MatGetRow_SeqBAIJ" 196 int MatGetRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 197 { 198 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 199 int itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*idx_i,bs2; 200 Scalar *aa,*v_i,*aa_i; 201 202 PetscFunctionBegin; 203 bs = a->bs; 204 ai = a->i; 205 aj = a->j; 206 aa = a->a; 207 bs2 = a->bs2; 208 209 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 210 211 bn = row/bs; /* Block number */ 212 bp = row % bs; /* Block Position */ 213 M = ai[bn+1] - ai[bn]; 214 *nz = bs*M; 215 216 if (v) { 217 *v = 0; 218 if (*nz) { 219 *v = (Scalar *) PetscMalloc( (*nz)*sizeof(Scalar) ); CHKPTRQ(*v); 220 for ( i=0; i<M; i++ ) { /* for each block in the block row */ 221 v_i = *v + i*bs; 222 aa_i = aa + bs2*(ai[bn] + i); 223 for ( j=bp,k=0; j<bs2; j+=bs,k++ ) {v_i[k] = aa_i[j];} 224 } 225 } 226 } 227 228 if (idx) { 229 *idx = 0; 230 if (*nz) { 231 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 232 for ( i=0; i<M; i++ ) { /* for each block in the block row */ 233 idx_i = *idx + i*bs; 234 itmp = bs*aj[ai[bn] + i]; 235 for ( j=0; j<bs; j++ ) {idx_i[j] = itmp++;} 236 } 237 } 238 } 239 PetscFunctionReturn(0); 240 } 241 242 #undef __FUNC__ 243 #define __FUNC__ "MatRestoreRow_SeqBAIJ" 244 int MatRestoreRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 245 { 246 PetscFunctionBegin; 247 if (idx) {if (*idx) PetscFree(*idx);} 248 if (v) {if (*v) PetscFree(*v);} 249 PetscFunctionReturn(0); 250 } 251 252 #undef __FUNC__ 253 #define __FUNC__ "MatTranspose_SeqBAIJ" 254 int MatTranspose_SeqBAIJ(Mat A,Mat *B) 255 { 256 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 257 Mat C; 258 int i,j,k,ierr,*aj=a->j,*ai=a->i,bs=a->bs,mbs=a->mbs,nbs=a->nbs,len,*col; 259 int *rows,*cols,bs2=a->bs2; 260 Scalar *array=a->a; 261 262 PetscFunctionBegin; 263 if (B==PETSC_NULL && mbs!=nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Square matrix only for in-place"); 264 col = (int *) PetscMalloc((1+nbs)*sizeof(int)); CHKPTRQ(col); 265 PetscMemzero(col,(1+nbs)*sizeof(int)); 266 267 for ( i=0; i<ai[mbs]; i++ ) col[aj[i]] += 1; 268 ierr = MatCreateSeqBAIJ(A->comm,bs,a->n,a->m,PETSC_NULL,col,&C); CHKERRQ(ierr); 269 PetscFree(col); 270 rows = (int *) PetscMalloc(2*bs*sizeof(int)); CHKPTRQ(rows); 271 cols = rows + bs; 272 for ( i=0; i<mbs; i++ ) { 273 cols[0] = i*bs; 274 for (k=1; k<bs; k++ ) cols[k] = cols[k-1] + 1; 275 len = ai[i+1] - ai[i]; 276 for ( j=0; j<len; j++ ) { 277 rows[0] = (*aj++)*bs; 278 for (k=1; k<bs; k++ ) rows[k] = rows[k-1] + 1; 279 ierr = MatSetValues(C,bs,rows,bs,cols,array,INSERT_VALUES); CHKERRQ(ierr); 280 array += bs2; 281 } 282 } 283 PetscFree(rows); 284 285 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 286 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 287 288 if (B != PETSC_NULL) { 289 *B = C; 290 } else { 291 PetscOps *Abops; 292 MatOps Aops; 293 294 /* This isn't really an in-place transpose */ 295 PetscFree(a->a); 296 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 297 if (a->diag) PetscFree(a->diag); 298 if (a->ilen) PetscFree(a->ilen); 299 if (a->imax) PetscFree(a->imax); 300 if (a->solve_work) PetscFree(a->solve_work); 301 PetscFree(a); 302 303 304 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 305 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 306 307 /* 308 This is horrible, horrible code. We need to keep the 309 A pointers for the bops and ops but copy everything 310 else from C. 311 */ 312 Abops = A->bops; 313 Aops = A->ops; 314 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 315 A->bops = Abops; 316 A->ops = Aops; 317 A->qlist = 0; 318 319 PetscHeaderDestroy(C); 320 } 321 PetscFunctionReturn(0); 322 } 323 324 325 326 327 #undef __FUNC__ 328 #define __FUNC__ "MatView_SeqBAIJ_Binary" 329 static int MatView_SeqBAIJ_Binary(Mat A,Viewer viewer) 330 { 331 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 332 int i, fd, *col_lens, ierr, bs = a->bs,count,*jj,j,k,l,bs2=a->bs2; 333 Scalar *aa; 334 335 PetscFunctionBegin; 336 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 337 col_lens = (int *) PetscMalloc((4+a->m)*sizeof(int));CHKPTRQ(col_lens); 338 col_lens[0] = MAT_COOKIE; 339 340 col_lens[1] = a->m; 341 col_lens[2] = a->n; 342 col_lens[3] = a->nz*bs2; 343 344 /* store lengths of each row and write (including header) to file */ 345 count = 0; 346 for ( i=0; i<a->mbs; i++ ) { 347 for ( j=0; j<bs; j++ ) { 348 col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]); 349 } 350 } 351 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,PETSC_INT,1); CHKERRQ(ierr); 352 PetscFree(col_lens); 353 354 /* store column indices (zero start index) */ 355 jj = (int *) PetscMalloc( (a->nz+1)*bs2*sizeof(int) ); CHKPTRQ(jj); 356 count = 0; 357 for ( i=0; i<a->mbs; i++ ) { 358 for ( j=0; j<bs; j++ ) { 359 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 360 for ( l=0; l<bs; l++ ) { 361 jj[count++] = bs*a->j[k] + l; 362 } 363 } 364 } 365 } 366 ierr = PetscBinaryWrite(fd,jj,bs2*a->nz,PETSC_INT,0); CHKERRQ(ierr); 367 PetscFree(jj); 368 369 /* store nonzero values */ 370 aa = (Scalar *) PetscMalloc((a->nz+1)*bs2*sizeof(Scalar)); CHKPTRQ(aa); 371 count = 0; 372 for ( i=0; i<a->mbs; i++ ) { 373 for ( j=0; j<bs; j++ ) { 374 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 375 for ( l=0; l<bs; l++ ) { 376 aa[count++] = a->a[bs2*k + l*bs + j]; 377 } 378 } 379 } 380 } 381 ierr = PetscBinaryWrite(fd,aa,bs2*a->nz,PETSC_SCALAR,0); CHKERRQ(ierr); 382 PetscFree(aa); 383 PetscFunctionReturn(0); 384 } 385 386 #undef __FUNC__ 387 #define __FUNC__ "MatView_SeqBAIJ_ASCII" 388 static int MatView_SeqBAIJ_ASCII(Mat A,Viewer viewer) 389 { 390 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 391 int ierr, i,j,format,bs = a->bs,k,l,bs2=a->bs2; 392 FILE *fd; 393 char *outputname; 394 395 PetscFunctionBegin; 396 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 397 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 398 ierr = ViewerGetFormat(viewer,&format); 399 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 400 fprintf(fd," block size is %d\n",bs); 401 } 402 else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 403 SETERRQ(PETSC_ERR_SUP,0,"Matlab format not supported"); 404 } 405 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 406 for ( i=0; i<a->mbs; i++ ) { 407 for ( j=0; j<bs; j++ ) { 408 fprintf(fd,"row %d:",i*bs+j); 409 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 410 for ( l=0; l<bs; l++ ) { 411 #if defined(USE_PETSC_COMPLEX) 412 if (PetscImaginary(a->a[bs2*k + l*bs + j]) > 0.0 && PetscReal(a->a[bs2*k + l*bs + j]) != 0.0) 413 fprintf(fd," %d %g + %g i",bs*a->j[k]+l, 414 PetscReal(a->a[bs2*k + l*bs + j]),PetscImaginary(a->a[bs2*k + l*bs + j])); 415 else if (PetscImaginary(a->a[bs2*k + l*bs + j]) < 0.0 && PetscReal(a->a[bs2*k + l*bs + j]) != 0.0) 416 fprintf(fd," %d %g - %g i",bs*a->j[k]+l, 417 PetscReal(a->a[bs2*k + l*bs + j]),-PetscImaginary(a->a[bs2*k + l*bs + j])); 418 else if (PetscReal(a->a[bs2*k + l*bs + j]) != 0.0) 419 fprintf(fd," %d %g ",bs*a->j[k]+l,PetscReal(a->a[bs2*k + l*bs + j])); 420 #else 421 if (a->a[bs2*k + l*bs + j] != 0.0) 422 fprintf(fd," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]); 423 #endif 424 } 425 } 426 fprintf(fd,"\n"); 427 } 428 } 429 } 430 else { 431 for ( i=0; i<a->mbs; i++ ) { 432 for ( j=0; j<bs; j++ ) { 433 fprintf(fd,"row %d:",i*bs+j); 434 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 435 for ( l=0; l<bs; l++ ) { 436 #if defined(USE_PETSC_COMPLEX) 437 if (PetscImaginary(a->a[bs2*k + l*bs + j]) > 0.0) { 438 fprintf(fd," %d %g + %g i",bs*a->j[k]+l, 439 PetscReal(a->a[bs2*k + l*bs + j]),PetscImaginary(a->a[bs2*k + l*bs + j])); 440 } 441 else if (PetscImaginary(a->a[bs2*k + l*bs + j]) < 0.0) { 442 fprintf(fd," %d %g - %g i",bs*a->j[k]+l, 443 PetscReal(a->a[bs2*k + l*bs + j]),-PetscImaginary(a->a[bs2*k + l*bs + j])); 444 } 445 else { 446 fprintf(fd," %d %g ",bs*a->j[k]+l,PetscReal(a->a[bs2*k + l*bs + j])); 447 } 448 #else 449 fprintf(fd," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]); 450 #endif 451 } 452 } 453 fprintf(fd,"\n"); 454 } 455 } 456 } 457 fflush(fd); 458 PetscFunctionReturn(0); 459 } 460 461 #undef __FUNC__ 462 #define __FUNC__ "MatView_SeqBAIJ_Draw_Zoom" 463 static int MatView_SeqBAIJ_Draw_Zoom(Draw draw,void *Aa) 464 { 465 Mat A = (Mat) Aa; 466 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *) A->data; 467 int row,ierr,i,j,k,l,mbs=a->mbs,color,bs=a->bs,bs2=a->bs2,rank; 468 double xl,yl,xr,yr,yc,x_l,x_r,y_l,y_r; 469 Scalar *aa; 470 PetscTruth isnull; 471 MPI_Comm comm; 472 Viewer viewer; 473 474 PetscFunctionBegin; 475 /* 476 This is nasty. If this is called from an originally parallel matrix 477 then all processes call this, but only the first has the matrix so the 478 rest should return immediately. 479 */ 480 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 481 MPI_Comm_rank(comm,&rank); 482 if (rank) PetscFunctionReturn(0); 483 484 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*) &viewer);CHKERRQ(ierr); 485 486 ierr = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr); CHKERRQ(ierr); 487 488 /* loop over matrix elements drawing boxes */ 489 color = DRAW_BLUE; 490 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 491 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 492 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 493 x_l = a->j[j]*bs; x_r = x_l + 1.0; 494 aa = a->a + j*bs2; 495 for ( k=0; k<bs; k++ ) { 496 for ( l=0; l<bs; l++ ) { 497 if (PetscReal(*aa++) >= 0.) continue; 498 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 499 } 500 } 501 } 502 } 503 color = DRAW_CYAN; 504 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 505 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 506 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 507 x_l = a->j[j]*bs; x_r = x_l + 1.0; 508 aa = a->a + j*bs2; 509 for ( k=0; k<bs; k++ ) { 510 for ( l=0; l<bs; l++ ) { 511 if (PetscReal(*aa++) != 0.) continue; 512 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 513 } 514 } 515 } 516 } 517 518 color = DRAW_RED; 519 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 520 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 521 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 522 x_l = a->j[j]*bs; x_r = x_l + 1.0; 523 aa = a->a + j*bs2; 524 for ( k=0; k<bs; k++ ) { 525 for ( l=0; l<bs; l++ ) { 526 if (PetscReal(*aa++) <= 0.) continue; 527 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 528 } 529 } 530 } 531 } 532 PetscFunctionReturn(0); 533 } 534 535 #undef __FUNC__ 536 #define __FUNC__ "MatView_SeqBAIJ_Draw" 537 static int MatView_SeqBAIJ_Draw(Mat A,Viewer viewer) 538 { 539 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *) A->data; 540 int ierr; 541 double xl,yl,xr,yr,w,h; 542 Draw draw; 543 PetscTruth isnull; 544 545 PetscFunctionBegin; 546 547 ierr = ViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 548 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 549 550 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 551 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 552 xr += w; yr += h; xl = -w; yl = -h; 553 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 554 ierr = DrawZoom(draw,MatView_SeqBAIJ_Draw_Zoom,A); CHKERRQ(ierr); 555 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 556 PetscFunctionReturn(0); 557 } 558 559 #undef __FUNC__ 560 #define __FUNC__ "MatView_SeqBAIJ" 561 int MatView_SeqBAIJ(Mat A,Viewer viewer) 562 { 563 ViewerType vtype; 564 int ierr; 565 566 PetscFunctionBegin; 567 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 568 if (!PetscStrcmp(vtype,MATLAB_VIEWER)) { 569 SETERRQ(PETSC_ERR_SUP,0,"Matlab viewer not supported"); 570 } else if (!PetscStrcmp(vtype,ASCII_VIEWER)){ 571 ierr = MatView_SeqBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 572 } else if (!PetscStrcmp(vtype,BINARY_VIEWER)) { 573 ierr = MatView_SeqBAIJ_Binary(A,viewer);CHKERRQ(ierr); 574 } else if (!PetscStrcmp(vtype,DRAW_VIEWER)) { 575 ierr = MatView_SeqBAIJ_Draw(A,viewer);CHKERRQ(ierr); 576 } else { 577 SETERRQ(1,1,"Viewer type not supported by PETSc object"); 578 } 579 PetscFunctionReturn(0); 580 } 581 582 583 #undef __FUNC__ 584 #define __FUNC__ "MatGetValues_SeqBAIJ" 585 int MatGetValues_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 586 { 587 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 588 int *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 589 int *ai = a->i, *ailen = a->ilen; 590 int brow,bcol,ridx,cidx,bs=a->bs,bs2=a->bs2; 591 Scalar *ap, *aa = a->a, zero = 0.0; 592 593 PetscFunctionBegin; 594 for ( k=0; k<m; k++ ) { /* loop over rows */ 595 row = im[k]; brow = row/bs; 596 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 597 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 598 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 599 nrow = ailen[brow]; 600 for ( l=0; l<n; l++ ) { /* loop over columns */ 601 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 602 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 603 col = in[l] ; 604 bcol = col/bs; 605 cidx = col%bs; 606 ridx = row%bs; 607 high = nrow; 608 low = 0; /* assume unsorted */ 609 while (high-low > 5) { 610 t = (low+high)/2; 611 if (rp[t] > bcol) high = t; 612 else low = t; 613 } 614 for ( i=low; i<high; i++ ) { 615 if (rp[i] > bcol) break; 616 if (rp[i] == bcol) { 617 *v++ = ap[bs2*i+bs*cidx+ridx]; 618 goto finished; 619 } 620 } 621 *v++ = zero; 622 finished:; 623 } 624 } 625 PetscFunctionReturn(0); 626 } 627 628 629 #undef __FUNC__ 630 #define __FUNC__ "MatSetValuesBlocked_SeqBAIJ" 631 int MatSetValuesBlocked_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 632 { 633 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 634 int *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 635 int *imax=a->imax,*ai=a->i,*ailen=a->ilen, roworiented=a->roworiented; 636 int *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=a->bs,stepval; 637 Scalar *ap,*value=v,*aa=a->a,*bap; 638 639 PetscFunctionBegin; 640 if (roworiented) { 641 stepval = (n-1)*bs; 642 } else { 643 stepval = (m-1)*bs; 644 } 645 for ( k=0; k<m; k++ ) { /* loop over added rows */ 646 row = im[k]; 647 #if defined(USE_PETSC_BOPT_g) 648 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 649 if (row >= a->mbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 650 #endif 651 rp = aj + ai[row]; 652 ap = aa + bs2*ai[row]; 653 rmax = imax[row]; 654 nrow = ailen[row]; 655 low = 0; 656 for ( l=0; l<n; l++ ) { /* loop over added columns */ 657 #if defined(USE_PETSC_BOPT_g) 658 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 659 if (in[l] >= a->nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 660 #endif 661 col = in[l]; 662 if (roworiented) { 663 value = v + k*(stepval+bs)*bs + l*bs; 664 } else { 665 value = v + l*(stepval+bs)*bs + k*bs; 666 } 667 if (!sorted) low = 0; high = nrow; 668 while (high-low > 7) { 669 t = (low+high)/2; 670 if (rp[t] > col) high = t; 671 else low = t; 672 } 673 for ( i=low; i<high; i++ ) { 674 if (rp[i] > col) break; 675 if (rp[i] == col) { 676 bap = ap + bs2*i; 677 if (roworiented) { 678 if (is == ADD_VALUES) { 679 for ( ii=0; ii<bs; ii++,value+=stepval ) { 680 for (jj=ii; jj<bs2; jj+=bs ) { 681 bap[jj] += *value++; 682 } 683 } 684 } else { 685 for ( ii=0; ii<bs; ii++,value+=stepval ) { 686 for (jj=ii; jj<bs2; jj+=bs ) { 687 bap[jj] = *value++; 688 } 689 } 690 } 691 } else { 692 if (is == ADD_VALUES) { 693 for ( ii=0; ii<bs; ii++,value+=stepval ) { 694 for (jj=0; jj<bs; jj++ ) { 695 *bap++ += *value++; 696 } 697 } 698 } else { 699 for ( ii=0; ii<bs; ii++,value+=stepval ) { 700 for (jj=0; jj<bs; jj++ ) { 701 *bap++ = *value++; 702 } 703 } 704 } 705 } 706 goto noinsert2; 707 } 708 } 709 if (nonew == 1) goto noinsert2; 710 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 711 if (nrow >= rmax) { 712 /* there is no extra room in row, therefore enlarge */ 713 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 714 Scalar *new_a; 715 716 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 717 718 /* malloc new storage space */ 719 len = new_nz*(sizeof(int)+bs2*sizeof(Scalar))+(a->mbs+1)*sizeof(int); 720 new_a = (Scalar *) PetscMalloc( len ); CHKPTRQ(new_a); 721 new_j = (int *) (new_a + bs2*new_nz); 722 new_i = new_j + new_nz; 723 724 /* copy over old data into new slots */ 725 for ( ii=0; ii<row+1; ii++ ) {new_i[ii] = ai[ii];} 726 for ( ii=row+1; ii<a->mbs+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 727 PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(int)); 728 len = (new_nz - CHUNKSIZE - ai[row] - nrow); 729 PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(int)); 730 PetscMemcpy(new_a,aa,(ai[row]+nrow)*bs2*sizeof(Scalar)); 731 PetscMemzero(new_a+bs2*(ai[row]+nrow),bs2*CHUNKSIZE*sizeof(Scalar)); 732 PetscMemcpy(new_a+bs2*(ai[row]+nrow+CHUNKSIZE),aa+bs2*(ai[row]+nrow),bs2*len*sizeof(Scalar)); 733 /* free up old matrix storage */ 734 PetscFree(a->a); 735 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 736 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 737 a->singlemalloc = 1; 738 739 rp = aj + ai[row]; ap = aa + bs2*ai[row]; 740 rmax = imax[row] = imax[row] + CHUNKSIZE; 741 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(Scalar))); 742 a->maxnz += bs2*CHUNKSIZE; 743 a->reallocs++; 744 a->nz++; 745 } 746 N = nrow++ - 1; 747 /* shift up all the later entries in this row */ 748 for ( ii=N; ii>=i; ii-- ) { 749 rp[ii+1] = rp[ii]; 750 PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(Scalar)); 751 } 752 if (N >= i) PetscMemzero(ap+bs2*i,bs2*sizeof(Scalar)); 753 rp[i] = col; 754 bap = ap + bs2*i; 755 if (roworiented) { 756 for ( ii=0; ii<bs; ii++,value+=stepval) { 757 for (jj=ii; jj<bs2; jj+=bs ) { 758 bap[jj] = *value++; 759 } 760 } 761 } else { 762 for ( ii=0; ii<bs; ii++,value+=stepval ) { 763 for (jj=0; jj<bs; jj++ ) { 764 *bap++ = *value++; 765 } 766 } 767 } 768 noinsert2:; 769 low = i; 770 } 771 ailen[row] = nrow; 772 } 773 PetscFunctionReturn(0); 774 } 775 776 777 #undef __FUNC__ 778 #define __FUNC__ "MatAssemblyEnd_SeqBAIJ" 779 int MatAssemblyEnd_SeqBAIJ(Mat A,MatAssemblyType mode) 780 { 781 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 782 int fshift = 0,i,j,*ai = a->i, *aj = a->j, *imax = a->imax; 783 int m = a->m,*ip, N, *ailen = a->ilen; 784 int mbs = a->mbs, bs2 = a->bs2,rmax = 0; 785 Scalar *aa = a->a, *ap; 786 787 PetscFunctionBegin; 788 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 789 790 if (m) rmax = ailen[0]; 791 for ( i=1; i<mbs; i++ ) { 792 /* move each row back by the amount of empty slots (fshift) before it*/ 793 fshift += imax[i-1] - ailen[i-1]; 794 rmax = PetscMax(rmax,ailen[i]); 795 if (fshift) { 796 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 797 N = ailen[i]; 798 for ( j=0; j<N; j++ ) { 799 ip[j-fshift] = ip[j]; 800 PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(Scalar)); 801 } 802 } 803 ai[i] = ai[i-1] + ailen[i-1]; 804 } 805 if (mbs) { 806 fshift += imax[mbs-1] - ailen[mbs-1]; 807 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 808 } 809 /* reset ilen and imax for each row */ 810 for ( i=0; i<mbs; i++ ) { 811 ailen[i] = imax[i] = ai[i+1] - ai[i]; 812 } 813 a->nz = ai[mbs]; 814 815 /* diagonals may have moved, so kill the diagonal pointers */ 816 if (fshift && a->diag) { 817 PetscFree(a->diag); 818 PLogObjectMemory(A,-(m+1)*sizeof(int)); 819 a->diag = 0; 820 } 821 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Matrix size: %d X %d, block size %d; storage space: %d unneeded, %d used\n", 822 m,a->n,a->bs,fshift*bs2,a->nz*bs2); 823 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Number of mallocs during MatSetValues is %d\n", 824 a->reallocs); 825 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Most nonzeros blocks in any row is %d\n",rmax); 826 a->reallocs = 0; 827 A->info.nz_unneeded = (double)fshift*bs2; 828 829 PetscFunctionReturn(0); 830 } 831 832 833 /* idx should be of length atlease bs */ 834 #undef __FUNC__ 835 #define __FUNC__ "MatZeroRows_SeqBAIJ_Check_Block" 836 static int MatZeroRows_SeqBAIJ_Check_Block(int *idx, int bs, PetscTruth *flg) 837 { 838 int i,row; 839 840 PetscFunctionBegin; 841 row = idx[0]; 842 if (row%bs!=0) { *flg = PETSC_FALSE; PetscFunctionReturn(0); } 843 844 for ( i=1; i<bs; i++ ) { 845 if (row+i != idx[i]) { *flg = PETSC_FALSE; PetscFunctionReturn(0); } 846 } 847 *flg = PETSC_TRUE; 848 PetscFunctionReturn(0); 849 } 850 851 #undef __FUNC__ 852 #define __FUNC__ "MatZeroRows_SeqBAIJ" 853 int MatZeroRows_SeqBAIJ(Mat A,IS is, Scalar *diag) 854 { 855 Mat_SeqBAIJ *baij=(Mat_SeqBAIJ*)A->data; 856 IS is_local; 857 int ierr,i,j,count,m=baij->m,is_n,*is_idx,*rows,bs=baij->bs,bs2=baij->bs2; 858 PetscTruth flg; 859 Scalar zero = 0.0,*aa; 860 861 PetscFunctionBegin; 862 /* Make a copy of the IS and sort it */ 863 ierr = ISGetSize(is,&is_n);CHKERRQ(ierr); 864 ierr = ISGetIndices(is,&is_idx);CHKERRQ(ierr); 865 ierr = ISCreateGeneral(A->comm,is_n,is_idx,&is_local); CHKERRQ(ierr); 866 ierr = ISSort(is_local); CHKERRQ(ierr); 867 ierr = ISGetIndices(is_local,&rows); CHKERRQ(ierr); 868 869 i = 0; 870 while (i < is_n) { 871 if (rows[i]<0 || rows[i]>m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 872 flg = PETSC_FALSE; 873 if (i+bs <= is_n) {ierr = MatZeroRows_SeqBAIJ_Check_Block(rows+i,bs,&flg); CHKERRQ(ierr); } 874 count = (baij->i[rows[i]/bs +1] - baij->i[rows[i]/bs])*bs; 875 aa = baij->a + baij->i[rows[i]/bs]*bs2 + (rows[i]%bs); 876 if (flg) { /* There exists a block of rows to be Zerowed */ 877 if (baij->ilen[rows[i]/bs] > 0) { 878 PetscMemzero(aa,count*bs*sizeof(Scalar)); 879 baij->ilen[rows[i]/bs] = 1; 880 baij->j[baij->i[rows[i]/bs]] = rows[i]/bs; 881 } 882 i += bs; 883 } else { /* Zero out only the requested row */ 884 for ( j=0; j<count; j++ ) { 885 aa[0] = zero; 886 aa+=bs; 887 } 888 i++; 889 } 890 } 891 if (diag) { 892 for ( j=0; j<is_n; j++ ) { 893 ierr = (*A->ops->setvalues)(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); 894 /* ierr = MatSetValues(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); */ 895 } 896 } 897 ierr = ISRestoreIndices(is,&is_idx); CHKERRQ(ierr); 898 ierr = ISRestoreIndices(is_local,&rows); CHKERRQ(ierr); 899 ierr = ISDestroy(is_local); CHKERRQ(ierr); 900 ierr = MatAssemblyEnd_SeqBAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 901 PetscFunctionReturn(0); 902 } 903 904 #undef __FUNC__ 905 #define __FUNC__ "MatSetValues_SeqBAIJ" 906 int MatSetValues_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 907 { 908 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 909 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 910 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 911 int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol; 912 int ridx,cidx,bs2=a->bs2; 913 Scalar *ap,value,*aa=a->a,*bap; 914 915 PetscFunctionBegin; 916 for ( k=0; k<m; k++ ) { /* loop over added rows */ 917 row = im[k]; brow = row/bs; 918 #if defined(USE_PETSC_BOPT_g) 919 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 920 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 921 #endif 922 rp = aj + ai[brow]; 923 ap = aa + bs2*ai[brow]; 924 rmax = imax[brow]; 925 nrow = ailen[brow]; 926 low = 0; 927 for ( l=0; l<n; l++ ) { /* loop over added columns */ 928 #if defined(USE_PETSC_BOPT_g) 929 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 930 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 931 #endif 932 col = in[l]; bcol = col/bs; 933 ridx = row % bs; cidx = col % bs; 934 if (roworiented) { 935 value = *v++; 936 } else { 937 value = v[k + l*m]; 938 } 939 if (!sorted) low = 0; high = nrow; 940 while (high-low > 7) { 941 t = (low+high)/2; 942 if (rp[t] > bcol) high = t; 943 else low = t; 944 } 945 for ( i=low; i<high; i++ ) { 946 if (rp[i] > bcol) break; 947 if (rp[i] == bcol) { 948 bap = ap + bs2*i + bs*cidx + ridx; 949 if (is == ADD_VALUES) *bap += value; 950 else *bap = value; 951 goto noinsert1; 952 } 953 } 954 if (nonew == 1) goto noinsert1; 955 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 956 if (nrow >= rmax) { 957 /* there is no extra room in row, therefore enlarge */ 958 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 959 Scalar *new_a; 960 961 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 962 963 /* Malloc new storage space */ 964 len = new_nz*(sizeof(int)+bs2*sizeof(Scalar))+(a->mbs+1)*sizeof(int); 965 new_a = (Scalar *) PetscMalloc( len ); CHKPTRQ(new_a); 966 new_j = (int *) (new_a + bs2*new_nz); 967 new_i = new_j + new_nz; 968 969 /* copy over old data into new slots */ 970 for ( ii=0; ii<brow+1; ii++ ) {new_i[ii] = ai[ii];} 971 for ( ii=brow+1; ii<a->mbs+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 972 PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int)); 973 len = (new_nz - CHUNKSIZE - ai[brow] - nrow); 974 PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow, 975 len*sizeof(int)); 976 PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(Scalar)); 977 PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(Scalar)); 978 PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE), 979 aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(Scalar)); 980 /* free up old matrix storage */ 981 PetscFree(a->a); 982 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 983 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 984 a->singlemalloc = 1; 985 986 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 987 rmax = imax[brow] = imax[brow] + CHUNKSIZE; 988 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(Scalar))); 989 a->maxnz += bs2*CHUNKSIZE; 990 a->reallocs++; 991 a->nz++; 992 } 993 N = nrow++ - 1; 994 /* shift up all the later entries in this row */ 995 for ( ii=N; ii>=i; ii-- ) { 996 rp[ii+1] = rp[ii]; 997 PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(Scalar)); 998 } 999 if (N>=i) PetscMemzero(ap+bs2*i,bs2*sizeof(Scalar)); 1000 rp[i] = bcol; 1001 ap[bs2*i + bs*cidx + ridx] = value; 1002 noinsert1:; 1003 low = i; 1004 } 1005 ailen[brow] = nrow; 1006 } 1007 PetscFunctionReturn(0); 1008 } 1009 1010 extern int MatLUFactorSymbolic_SeqBAIJ(Mat,IS,IS,double,Mat*); 1011 extern int MatLUFactor_SeqBAIJ(Mat,IS,IS,double); 1012 extern int MatIncreaseOverlap_SeqBAIJ(Mat,int,IS*,int); 1013 extern int MatGetSubMatrix_SeqBAIJ(Mat,IS,IS,int,MatGetSubMatrixCall,Mat*); 1014 extern int MatGetSubMatrices_SeqBAIJ(Mat,int,IS*,IS*,MatGetSubMatrixCall,Mat**); 1015 extern int MatMultTrans_SeqBAIJ(Mat,Vec,Vec); 1016 extern int MatMultTransAdd_SeqBAIJ(Mat,Vec,Vec,Vec); 1017 extern int MatScale_SeqBAIJ(Scalar*,Mat); 1018 extern int MatNorm_SeqBAIJ(Mat,NormType,double *); 1019 extern int MatEqual_SeqBAIJ(Mat,Mat, PetscTruth*); 1020 extern int MatGetDiagonal_SeqBAIJ(Mat,Vec); 1021 extern int MatDiagonalScale_SeqBAIJ(Mat,Vec,Vec); 1022 extern int MatGetInfo_SeqBAIJ(Mat,MatInfoType,MatInfo *); 1023 extern int MatZeroEntries_SeqBAIJ(Mat); 1024 1025 extern int MatSolve_SeqBAIJ_N(Mat,Vec,Vec); 1026 extern int MatSolve_SeqBAIJ_1(Mat,Vec,Vec); 1027 extern int MatSolve_SeqBAIJ_2(Mat,Vec,Vec); 1028 extern int MatSolve_SeqBAIJ_3(Mat,Vec,Vec); 1029 extern int MatSolve_SeqBAIJ_4(Mat,Vec,Vec); 1030 extern int MatSolve_SeqBAIJ_5(Mat,Vec,Vec); 1031 extern int MatSolve_SeqBAIJ_7(Mat,Vec,Vec); 1032 1033 extern int MatLUFactorNumeric_SeqBAIJ_N(Mat,Mat*); 1034 extern int MatLUFactorNumeric_SeqBAIJ_1(Mat,Mat*); 1035 extern int MatLUFactorNumeric_SeqBAIJ_2(Mat,Mat*); 1036 extern int MatLUFactorNumeric_SeqBAIJ_3(Mat,Mat*); 1037 extern int MatLUFactorNumeric_SeqBAIJ_4(Mat,Mat*); 1038 extern int MatLUFactorNumeric_SeqBAIJ_5(Mat,Mat*); 1039 1040 extern int MatMult_SeqBAIJ_1(Mat,Vec,Vec); 1041 extern int MatMult_SeqBAIJ_2(Mat,Vec,Vec); 1042 extern int MatMult_SeqBAIJ_3(Mat,Vec,Vec); 1043 extern int MatMult_SeqBAIJ_4(Mat,Vec,Vec); 1044 extern int MatMult_SeqBAIJ_5(Mat,Vec,Vec); 1045 extern int MatMult_SeqBAIJ_7(Mat,Vec,Vec); 1046 extern int MatMult_SeqBAIJ_N(Mat,Vec,Vec); 1047 1048 extern int MatMultAdd_SeqBAIJ_1(Mat,Vec,Vec,Vec); 1049 extern int MatMultAdd_SeqBAIJ_2(Mat,Vec,Vec,Vec); 1050 extern int MatMultAdd_SeqBAIJ_3(Mat,Vec,Vec,Vec); 1051 extern int MatMultAdd_SeqBAIJ_4(Mat,Vec,Vec,Vec); 1052 extern int MatMultAdd_SeqBAIJ_5(Mat,Vec,Vec,Vec); 1053 extern int MatMultAdd_SeqBAIJ_7(Mat,Vec,Vec,Vec); 1054 extern int MatMultAdd_SeqBAIJ_N(Mat,Vec,Vec,Vec); 1055 1056 #undef __FUNC__ 1057 #define __FUNC__ "MatILUFactor_SeqBAIJ" 1058 int MatILUFactor_SeqBAIJ(Mat inA,IS row,IS col,double efill,int fill) 1059 { 1060 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) inA->data; 1061 Mat outA; 1062 int ierr; 1063 PetscTruth row_identity, col_identity; 1064 1065 PetscFunctionBegin; 1066 if (fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1067 ierr = ISIdentity(row,&row_identity); CHKERRQ(ierr); 1068 ierr = ISIdentity(col,&col_identity); CHKERRQ(ierr); 1069 if (!row_identity || !col_identity) { 1070 SETERRQ(1,1,"Row and column permutations must be identity"); 1071 } 1072 1073 outA = inA; 1074 inA->factor = FACTOR_LU; 1075 a->row = row; 1076 a->col = col; 1077 1078 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1079 ierr = ISInvertPermutation(col,&(a->icol)); CHKERRQ(ierr); 1080 PLogObjectParent(inA,a->icol); 1081 1082 if (!a->solve_work) { 1083 a->solve_work = (Scalar *) PetscMalloc((a->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1084 PLogObjectMemory(inA,(a->m+a->bs)*sizeof(Scalar)); 1085 } 1086 1087 if (!a->diag) { 1088 ierr = MatMarkDiag_SeqBAIJ(inA); CHKERRQ(ierr); 1089 } 1090 /* 1091 Blocksize 4 and 5 have a special faster factorization/solver for ILU(0) factorization 1092 with natural ordering 1093 */ 1094 if (a->bs == 4) { 1095 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 1096 inA->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 1097 PLogInfo(inA,"Using special in-place natural ordering factor and solve BS=4\n"); 1098 } else if (a->bs == 5) { 1099 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering; 1100 inA->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering; 1101 PLogInfo(inA,"Using special in-place natural ordering factor and solve BS=5\n"); 1102 } 1103 1104 ierr = MatLUFactorNumeric(inA,&outA); CHKERRQ(ierr); 1105 1106 1107 PetscFunctionReturn(0); 1108 } 1109 #undef __FUNC__ 1110 #define __FUNC__ "MatPrintHelp_SeqBAIJ" 1111 int MatPrintHelp_SeqBAIJ(Mat A) 1112 { 1113 static int called = 0; 1114 MPI_Comm comm = A->comm; 1115 1116 PetscFunctionBegin; 1117 if (called) {PetscFunctionReturn(0);} else called = 1; 1118 (*PetscHelpPrintf)(comm," Options for MATSEQBAIJ and MATMPIBAIJ matrix formats (the defaults):\n"); 1119 (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n"); 1120 PetscFunctionReturn(0); 1121 } 1122 1123 EXTERN_C_BEGIN 1124 #undef __FUNC__ 1125 #define __FUNC__ "MatSeqBAIJSetColumnIndices_SeqBAIJ" 1126 int MatSeqBAIJSetColumnIndices_SeqBAIJ(Mat mat,int *indices) 1127 { 1128 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data; 1129 int i,nz,n; 1130 1131 PetscFunctionBegin; 1132 nz = baij->maxnz; 1133 n = baij->n; 1134 for (i=0; i<nz; i++) { 1135 baij->j[i] = indices[i]; 1136 } 1137 baij->nz = nz; 1138 for ( i=0; i<n; i++ ) { 1139 baij->ilen[i] = baij->imax[i]; 1140 } 1141 1142 PetscFunctionReturn(0); 1143 } 1144 EXTERN_C_END 1145 1146 #undef __FUNC__ 1147 #define __FUNC__ "MatSeqBAIJSetColumnIndices" 1148 /*@ 1149 MatSeqBAIJSetColumnIndices - Set the column indices for all the rows 1150 in the matrix. 1151 1152 Input Parameters: 1153 + mat - the SeqBAIJ matrix 1154 - indices - the column indices 1155 1156 Notes: 1157 This can be called if you have precomputed the nonzero structure of the 1158 matrix and want to provide it to the matrix object to improve the performance 1159 of the MatSetValues() operation. 1160 1161 You MUST have set the correct numbers of nonzeros per row in the call to 1162 MatCreateSeqBAIJ(). 1163 1164 MUST be called before any calls to MatSetValues(); 1165 1166 @*/ 1167 int MatSeqBAIJSetColumnIndices(Mat mat,int *indices) 1168 { 1169 int ierr,(*f)(Mat,int *); 1170 1171 PetscFunctionBegin; 1172 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1173 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqBAIJSetColumnIndices_C",(void **)&f); 1174 CHKERRQ(ierr); 1175 if (f) { 1176 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1177 } else { 1178 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1179 } 1180 PetscFunctionReturn(0); 1181 } 1182 1183 /* -------------------------------------------------------------------*/ 1184 static struct _MatOps MatOps_Values = {MatSetValues_SeqBAIJ, 1185 MatGetRow_SeqBAIJ, 1186 MatRestoreRow_SeqBAIJ, 1187 MatMult_SeqBAIJ_N, 1188 MatMultAdd_SeqBAIJ_N, 1189 MatMultTrans_SeqBAIJ, 1190 MatMultTransAdd_SeqBAIJ, 1191 MatSolve_SeqBAIJ_N, 1192 0, 1193 0, 1194 0, 1195 MatLUFactor_SeqBAIJ, 1196 0, 1197 0, 1198 MatTranspose_SeqBAIJ, 1199 MatGetInfo_SeqBAIJ, 1200 MatEqual_SeqBAIJ, 1201 MatGetDiagonal_SeqBAIJ, 1202 MatDiagonalScale_SeqBAIJ, 1203 MatNorm_SeqBAIJ, 1204 0, 1205 MatAssemblyEnd_SeqBAIJ, 1206 0, 1207 MatSetOption_SeqBAIJ, 1208 MatZeroEntries_SeqBAIJ, 1209 MatZeroRows_SeqBAIJ, 1210 MatLUFactorSymbolic_SeqBAIJ, 1211 MatLUFactorNumeric_SeqBAIJ_N, 1212 0, 1213 0, 1214 MatGetSize_SeqBAIJ, 1215 MatGetSize_SeqBAIJ, 1216 MatGetOwnershipRange_SeqBAIJ, 1217 MatILUFactorSymbolic_SeqBAIJ, 1218 0, 1219 0, 1220 0, 1221 MatDuplicate_SeqBAIJ, 1222 0, 1223 0, 1224 MatILUFactor_SeqBAIJ, 1225 0, 1226 0, 1227 MatGetSubMatrices_SeqBAIJ, 1228 MatIncreaseOverlap_SeqBAIJ, 1229 MatGetValues_SeqBAIJ, 1230 0, 1231 MatPrintHelp_SeqBAIJ, 1232 MatScale_SeqBAIJ, 1233 0, 1234 0, 1235 0, 1236 MatGetBlockSize_SeqBAIJ, 1237 MatGetRowIJ_SeqBAIJ, 1238 MatRestoreRowIJ_SeqBAIJ, 1239 0, 1240 0, 1241 0, 1242 0, 1243 0, 1244 0, 1245 MatSetValuesBlocked_SeqBAIJ, 1246 MatGetSubMatrix_SeqBAIJ, 1247 0, 1248 0, 1249 MatGetMaps_Petsc}; 1250 1251 #undef __FUNC__ 1252 #define __FUNC__ "MatCreateSeqBAIJ" 1253 /*@C 1254 MatCreateSeqBAIJ - Creates a sparse matrix in block AIJ (block 1255 compressed row) format. For good matrix assembly performance the 1256 user should preallocate the matrix storage by setting the parameter nz 1257 (or the array nzz). By setting these parameters accurately, performance 1258 during matrix assembly can be increased by more than a factor of 50. 1259 1260 Collective on MPI_Comm 1261 1262 Input Parameters: 1263 + comm - MPI communicator, set to PETSC_COMM_SELF 1264 . bs - size of block 1265 . m - number of rows 1266 . n - number of columns 1267 . nz - number of block nonzeros per block row (same for all rows) 1268 - nzz - array containing the number of block nonzeros in the various block rows 1269 (possibly different for each block row) or PETSC_NULL 1270 1271 Output Parameter: 1272 . A - the matrix 1273 1274 Options Database Keys: 1275 . -mat_no_unroll - uses code that does not unroll the loops in the 1276 block calculations (much slower) 1277 . -mat_block_size - size of the blocks to use 1278 1279 Notes: 1280 The block AIJ format is fully compatible with standard Fortran 77 1281 storage. That is, the stored row and column indices can begin at 1282 either one (as in Fortran) or zero. See the users' manual for details. 1283 1284 Specify the preallocated storage with either nz or nnz (not both). 1285 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1286 allocation. For additional details, see the users manual chapter on 1287 matrices. 1288 1289 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1290 @*/ 1291 int MatCreateSeqBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz, Mat *A) 1292 { 1293 Mat B; 1294 Mat_SeqBAIJ *b; 1295 int i,len,ierr,flg,mbs=m/bs,nbs=n/bs,bs2=bs*bs,size; 1296 1297 PetscFunctionBegin; 1298 MPI_Comm_size(comm,&size); 1299 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1300 1301 ierr = OptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,&flg);CHKERRQ(ierr); 1302 1303 if (mbs*bs!=m || nbs*bs!=n) { 1304 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Number rows, cols must be divisible by blocksize"); 1305 } 1306 1307 *A = 0; 1308 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQBAIJ,comm,MatDestroy,MatView); 1309 PLogObjectCreate(B); 1310 B->data = (void *) (b = PetscNew(Mat_SeqBAIJ)); CHKPTRQ(b); 1311 PetscMemzero(b,sizeof(Mat_SeqBAIJ)); 1312 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 1313 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg); CHKERRQ(ierr); 1314 if (!flg) { 1315 switch (bs) { 1316 case 1: 1317 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 1318 B->ops->solve = MatSolve_SeqBAIJ_1; 1319 B->ops->mult = MatMult_SeqBAIJ_1; 1320 B->ops->multadd = MatMultAdd_SeqBAIJ_1; 1321 break; 1322 case 2: 1323 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2; 1324 B->ops->solve = MatSolve_SeqBAIJ_2; 1325 B->ops->mult = MatMult_SeqBAIJ_2; 1326 B->ops->multadd = MatMultAdd_SeqBAIJ_2; 1327 break; 1328 case 3: 1329 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3; 1330 B->ops->solve = MatSolve_SeqBAIJ_3; 1331 B->ops->mult = MatMult_SeqBAIJ_3; 1332 B->ops->multadd = MatMultAdd_SeqBAIJ_3; 1333 break; 1334 case 4: 1335 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4; 1336 B->ops->solve = MatSolve_SeqBAIJ_4; 1337 B->ops->mult = MatMult_SeqBAIJ_4; 1338 B->ops->multadd = MatMultAdd_SeqBAIJ_4; 1339 break; 1340 case 5: 1341 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5; 1342 B->ops->solve = MatSolve_SeqBAIJ_5; 1343 B->ops->mult = MatMult_SeqBAIJ_5; 1344 B->ops->multadd = MatMultAdd_SeqBAIJ_5; 1345 break; 1346 case 7: 1347 B->ops->mult = MatMult_SeqBAIJ_7; 1348 B->ops->solve = MatSolve_SeqBAIJ_7; 1349 B->ops->multadd = MatMultAdd_SeqBAIJ_7; 1350 break; 1351 } 1352 } 1353 B->ops->destroy = MatDestroy_SeqBAIJ; 1354 B->ops->view = MatView_SeqBAIJ; 1355 B->factor = 0; 1356 B->lupivotthreshold = 1.0; 1357 B->mapping = 0; 1358 b->row = 0; 1359 b->col = 0; 1360 b->icol = 0; 1361 b->reallocs = 0; 1362 1363 b->m = m; B->m = m; B->M = m; 1364 b->n = n; B->n = n; B->N = n; 1365 1366 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1367 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 1368 1369 b->mbs = mbs; 1370 b->nbs = nbs; 1371 b->imax = (int *) PetscMalloc( (mbs+1)*sizeof(int) ); CHKPTRQ(b->imax); 1372 if (nnz == PETSC_NULL) { 1373 if (nz == PETSC_DEFAULT) nz = 5; 1374 else if (nz <= 0) nz = 1; 1375 for ( i=0; i<mbs; i++ ) b->imax[i] = nz; 1376 nz = nz*mbs; 1377 } else { 1378 nz = 0; 1379 for ( i=0; i<mbs; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1380 } 1381 1382 /* allocate the matrix space */ 1383 len = nz*sizeof(int) + nz*bs2*sizeof(Scalar) + (b->m+1)*sizeof(int); 1384 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1385 PetscMemzero(b->a,nz*bs2*sizeof(Scalar)); 1386 b->j = (int *) (b->a + nz*bs2); 1387 PetscMemzero(b->j,nz*sizeof(int)); 1388 b->i = b->j + nz; 1389 b->singlemalloc = 1; 1390 1391 b->i[0] = 0; 1392 for (i=1; i<mbs+1; i++) { 1393 b->i[i] = b->i[i-1] + b->imax[i-1]; 1394 } 1395 1396 /* b->ilen will count nonzeros in each block row so far. */ 1397 b->ilen = (int *) PetscMalloc((mbs+1)*sizeof(int)); 1398 PLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1399 for ( i=0; i<mbs; i++ ) { b->ilen[i] = 0;} 1400 1401 b->bs = bs; 1402 b->bs2 = bs2; 1403 b->mbs = mbs; 1404 b->nz = 0; 1405 b->maxnz = nz*bs2; 1406 b->sorted = 0; 1407 b->roworiented = 1; 1408 b->nonew = 0; 1409 b->diag = 0; 1410 b->solve_work = 0; 1411 b->mult_work = 0; 1412 b->spptr = 0; 1413 B->info.nz_unneeded = (double)b->maxnz; 1414 1415 *A = B; 1416 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1417 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1418 1419 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqBAIJSetColumnIndices_C", 1420 "MatSeqBAIJSetColumnIndices_SeqBAIJ", 1421 (void*)MatSeqBAIJSetColumnIndices_SeqBAIJ);CHKERRQ(ierr); 1422 1423 PetscFunctionReturn(0); 1424 } 1425 1426 #undef __FUNC__ 1427 #define __FUNC__ "MatDuplicate_SeqBAIJ" 1428 int MatDuplicate_SeqBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1429 { 1430 Mat C; 1431 Mat_SeqBAIJ *c,*a = (Mat_SeqBAIJ *) A->data; 1432 int i,len, mbs = a->mbs,nz = a->nz,bs2 =a->bs2,ierr; 1433 1434 PetscFunctionBegin; 1435 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,0,"Corrupt matrix"); 1436 1437 *B = 0; 1438 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQBAIJ,A->comm,MatDestroy,MatView); 1439 PLogObjectCreate(C); 1440 C->data = (void *) (c = PetscNew(Mat_SeqBAIJ)); CHKPTRQ(c); 1441 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 1442 C->ops->destroy = MatDestroy_SeqBAIJ; 1443 C->ops->view = MatView_SeqBAIJ; 1444 C->factor = A->factor; 1445 c->row = 0; 1446 c->col = 0; 1447 C->assembled = PETSC_TRUE; 1448 1449 c->m = C->m = a->m; 1450 c->n = C->n = a->n; 1451 C->M = a->m; 1452 C->N = a->n; 1453 1454 c->bs = a->bs; 1455 c->bs2 = a->bs2; 1456 c->mbs = a->mbs; 1457 c->nbs = a->nbs; 1458 1459 c->imax = (int *) PetscMalloc((mbs+1)*sizeof(int)); CHKPTRQ(c->imax); 1460 c->ilen = (int *) PetscMalloc((mbs+1)*sizeof(int)); CHKPTRQ(c->ilen); 1461 for ( i=0; i<mbs; i++ ) { 1462 c->imax[i] = a->imax[i]; 1463 c->ilen[i] = a->ilen[i]; 1464 } 1465 1466 /* allocate the matrix space */ 1467 c->singlemalloc = 1; 1468 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(Scalar) + sizeof(int)); 1469 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1470 c->j = (int *) (c->a + nz*bs2); 1471 c->i = c->j + nz; 1472 PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int)); 1473 if (mbs > 0) { 1474 PetscMemcpy(c->j,a->j,nz*sizeof(int)); 1475 if (cpvalues == MAT_COPY_VALUES) { 1476 PetscMemcpy(c->a,a->a,bs2*nz*sizeof(Scalar)); 1477 } else { 1478 PetscMemzero(c->a,bs2*nz*sizeof(Scalar)); 1479 } 1480 } 1481 1482 PLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1483 c->sorted = a->sorted; 1484 c->roworiented = a->roworiented; 1485 c->nonew = a->nonew; 1486 1487 if (a->diag) { 1488 c->diag = (int *) PetscMalloc( (mbs+1)*sizeof(int) ); CHKPTRQ(c->diag); 1489 PLogObjectMemory(C,(mbs+1)*sizeof(int)); 1490 for ( i=0; i<mbs; i++ ) { 1491 c->diag[i] = a->diag[i]; 1492 } 1493 } 1494 else c->diag = 0; 1495 c->nz = a->nz; 1496 c->maxnz = a->maxnz; 1497 c->solve_work = 0; 1498 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1499 c->mult_work = 0; 1500 *B = C; 1501 ierr = PetscObjectComposeFunction((PetscObject)C,"MatSeqBAIJSetColumnIndices_C", 1502 "MatSeqBAIJSetColumnIndices_SeqBAIJ", 1503 (void*)MatSeqBAIJSetColumnIndices_SeqBAIJ);CHKERRQ(ierr); 1504 PetscFunctionReturn(0); 1505 } 1506 1507 #undef __FUNC__ 1508 #define __FUNC__ "MatLoad_SeqBAIJ" 1509 int MatLoad_SeqBAIJ(Viewer viewer,MatType type,Mat *A) 1510 { 1511 Mat_SeqBAIJ *a; 1512 Mat B; 1513 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1,flg; 1514 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount; 1515 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1516 int *masked, nmask,tmp,bs2,ishift; 1517 Scalar *aa; 1518 MPI_Comm comm = ((PetscObject) viewer)->comm; 1519 1520 PetscFunctionBegin; 1521 ierr = OptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,&flg);CHKERRQ(ierr); 1522 bs2 = bs*bs; 1523 1524 MPI_Comm_size(comm,&size); 1525 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 1526 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1527 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 1528 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not Mat object"); 1529 M = header[1]; N = header[2]; nz = header[3]; 1530 1531 if (header[3] < 0) { 1532 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format, cannot load as SeqBAIJ"); 1533 } 1534 1535 if (M != N) SETERRQ(PETSC_ERR_SUP,0,"Can only do square matrices"); 1536 1537 /* 1538 This code adds extra rows to make sure the number of rows is 1539 divisible by the blocksize 1540 */ 1541 mbs = M/bs; 1542 extra_rows = bs - M + bs*(mbs); 1543 if (extra_rows == bs) extra_rows = 0; 1544 else mbs++; 1545 if (extra_rows) { 1546 PLogInfo(0,"MatLoad_SeqBAIJ:Padding loaded matrix to match blocksize\n"); 1547 } 1548 1549 /* read in row lengths */ 1550 rowlengths = (int*) PetscMalloc((M+extra_rows)*sizeof(int));CHKPTRQ(rowlengths); 1551 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 1552 for ( i=0; i<extra_rows; i++ ) rowlengths[M+i] = 1; 1553 1554 /* read in column indices */ 1555 jj = (int*) PetscMalloc( (nz+extra_rows)*sizeof(int) ); CHKPTRQ(jj); 1556 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT); CHKERRQ(ierr); 1557 for ( i=0; i<extra_rows; i++ ) jj[nz+i] = M+i; 1558 1559 /* loop over row lengths determining block row lengths */ 1560 browlengths = (int *) PetscMalloc(mbs*sizeof(int));CHKPTRQ(browlengths); 1561 PetscMemzero(browlengths,mbs*sizeof(int)); 1562 mask = (int *) PetscMalloc( 2*mbs*sizeof(int) ); CHKPTRQ(mask); 1563 PetscMemzero(mask,mbs*sizeof(int)); 1564 masked = mask + mbs; 1565 rowcount = 0; nzcount = 0; 1566 for ( i=0; i<mbs; i++ ) { 1567 nmask = 0; 1568 for ( j=0; j<bs; j++ ) { 1569 kmax = rowlengths[rowcount]; 1570 for ( k=0; k<kmax; k++ ) { 1571 tmp = jj[nzcount++]/bs; 1572 if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;} 1573 } 1574 rowcount++; 1575 } 1576 browlengths[i] += nmask; 1577 /* zero out the mask elements we set */ 1578 for ( j=0; j<nmask; j++ ) mask[masked[j]] = 0; 1579 } 1580 1581 /* create our matrix */ 1582 ierr = MatCreateSeqBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1583 B = *A; 1584 a = (Mat_SeqBAIJ *) B->data; 1585 1586 /* set matrix "i" values */ 1587 a->i[0] = 0; 1588 for ( i=1; i<= mbs; i++ ) { 1589 a->i[i] = a->i[i-1] + browlengths[i-1]; 1590 a->ilen[i-1] = browlengths[i-1]; 1591 } 1592 a->nz = 0; 1593 for ( i=0; i<mbs; i++ ) a->nz += browlengths[i]; 1594 1595 /* read in nonzero values */ 1596 aa = (Scalar *) PetscMalloc((nz+extra_rows)*sizeof(Scalar));CHKPTRQ(aa); 1597 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR); CHKERRQ(ierr); 1598 for ( i=0; i<extra_rows; i++ ) aa[nz+i] = 1.0; 1599 1600 /* set "a" and "j" values into matrix */ 1601 nzcount = 0; jcount = 0; 1602 for ( i=0; i<mbs; i++ ) { 1603 nzcountb = nzcount; 1604 nmask = 0; 1605 for ( j=0; j<bs; j++ ) { 1606 kmax = rowlengths[i*bs+j]; 1607 for ( k=0; k<kmax; k++ ) { 1608 tmp = jj[nzcount++]/bs; 1609 if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;} 1610 } 1611 rowcount++; 1612 } 1613 /* sort the masked values */ 1614 PetscSortInt(nmask,masked); 1615 1616 /* set "j" values into matrix */ 1617 maskcount = 1; 1618 for ( j=0; j<nmask; j++ ) { 1619 a->j[jcount++] = masked[j]; 1620 mask[masked[j]] = maskcount++; 1621 } 1622 /* set "a" values into matrix */ 1623 ishift = bs2*a->i[i]; 1624 for ( j=0; j<bs; j++ ) { 1625 kmax = rowlengths[i*bs+j]; 1626 for ( k=0; k<kmax; k++ ) { 1627 tmp = jj[nzcountb]/bs ; 1628 block = mask[tmp] - 1; 1629 point = jj[nzcountb] - bs*tmp; 1630 idx = ishift + bs2*block + j + bs*point; 1631 a->a[idx] = aa[nzcountb++]; 1632 } 1633 } 1634 /* zero out the mask elements we set */ 1635 for ( j=0; j<nmask; j++ ) mask[masked[j]] = 0; 1636 } 1637 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Bad binary matrix"); 1638 1639 PetscFree(rowlengths); 1640 PetscFree(browlengths); 1641 PetscFree(aa); 1642 PetscFree(jj); 1643 PetscFree(mask); 1644 1645 B->assembled = PETSC_TRUE; 1646 1647 ierr = MatView_Private(B); CHKERRQ(ierr); 1648 PetscFunctionReturn(0); 1649 } 1650 1651 1652 1653