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