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