1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: baij.c,v 1.123 1998/01/14 02:41:43 bsmith Exp bsmith $"; 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,"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 if (baij->ilen[rows[i]/bs] > 0) { 874 PetscMemzero(aa,count*bs*sizeof(Scalar)); 875 baij->ilen[rows[i]/bs] = 1; 876 baij->j[baij->i[rows[i]/bs]] = rows[i]/bs; 877 } 878 i += bs; 879 } else { /* Zero out only the requested row */ 880 for ( j=0; j<count; j++ ) { 881 aa[0] = zero; 882 aa+=bs; 883 } 884 i++; 885 } 886 } 887 if (diag) { 888 for ( j=0; j<is_n; j++ ) { 889 ierr = (*A->ops.setvalues)(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); 890 /* ierr = MatSetValues(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); */ 891 } 892 } 893 ierr = ISRestoreIndices(is,&is_idx); CHKERRQ(ierr); 894 ierr = ISRestoreIndices(is_local,&rows); CHKERRQ(ierr); 895 ierr = ISDestroy(is_local); CHKERRQ(ierr); 896 ierr = MatAssemblyEnd_SeqBAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 897 PetscFunctionReturn(0); 898 } 899 900 #undef __FUNC__ 901 #define __FUNC__ "MatSetValues_SeqBAIJ" 902 int MatSetValues_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 903 { 904 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 905 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 906 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 907 int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol; 908 int ridx,cidx,bs2=a->bs2; 909 Scalar *ap,value,*aa=a->a,*bap; 910 911 PetscFunctionBegin; 912 for ( k=0; k<m; k++ ) { /* loop over added rows */ 913 row = im[k]; brow = row/bs; 914 #if defined(USE_PETSC_BOPT_g) 915 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 916 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 917 #endif 918 rp = aj + ai[brow]; 919 ap = aa + bs2*ai[brow]; 920 rmax = imax[brow]; 921 nrow = ailen[brow]; 922 low = 0; 923 for ( l=0; l<n; l++ ) { /* loop over added columns */ 924 #if defined(USE_PETSC_BOPT_g) 925 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 926 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 927 #endif 928 col = in[l]; bcol = col/bs; 929 ridx = row % bs; cidx = col % bs; 930 if (roworiented) { 931 value = *v++; 932 } else { 933 value = v[k + l*m]; 934 } 935 if (!sorted) low = 0; high = nrow; 936 while (high-low > 7) { 937 t = (low+high)/2; 938 if (rp[t] > bcol) high = t; 939 else low = t; 940 } 941 for ( i=low; i<high; i++ ) { 942 if (rp[i] > bcol) break; 943 if (rp[i] == bcol) { 944 bap = ap + bs2*i + bs*cidx + ridx; 945 if (is == ADD_VALUES) *bap += value; 946 else *bap = value; 947 goto noinsert1; 948 } 949 } 950 if (nonew == 1) goto noinsert1; 951 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 952 if (nrow >= rmax) { 953 /* there is no extra room in row, therefore enlarge */ 954 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 955 Scalar *new_a; 956 957 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 958 959 /* Malloc new storage space */ 960 len = new_nz*(sizeof(int)+bs2*sizeof(Scalar))+(a->mbs+1)*sizeof(int); 961 new_a = (Scalar *) PetscMalloc( len ); CHKPTRQ(new_a); 962 new_j = (int *) (new_a + bs2*new_nz); 963 new_i = new_j + new_nz; 964 965 /* copy over old data into new slots */ 966 for ( ii=0; ii<brow+1; ii++ ) {new_i[ii] = ai[ii];} 967 for ( ii=brow+1; ii<a->mbs+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 968 PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int)); 969 len = (new_nz - CHUNKSIZE - ai[brow] - nrow); 970 PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow, 971 len*sizeof(int)); 972 PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(Scalar)); 973 PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(Scalar)); 974 PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE), 975 aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(Scalar)); 976 /* free up old matrix storage */ 977 PetscFree(a->a); 978 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 979 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 980 a->singlemalloc = 1; 981 982 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 983 rmax = imax[brow] = imax[brow] + CHUNKSIZE; 984 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(Scalar))); 985 a->maxnz += bs2*CHUNKSIZE; 986 a->reallocs++; 987 a->nz++; 988 } 989 N = nrow++ - 1; 990 /* shift up all the later entries in this row */ 991 for ( ii=N; ii>=i; ii-- ) { 992 rp[ii+1] = rp[ii]; 993 PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(Scalar)); 994 } 995 if (N>=i) PetscMemzero(ap+bs2*i,bs2*sizeof(Scalar)); 996 rp[i] = bcol; 997 ap[bs2*i + bs*cidx + ridx] = value; 998 noinsert1:; 999 low = i; 1000 } 1001 ailen[brow] = nrow; 1002 } 1003 PetscFunctionReturn(0); 1004 } 1005 1006 extern int MatLUFactorSymbolic_SeqBAIJ(Mat,IS,IS,double,Mat*); 1007 extern int MatLUFactor_SeqBAIJ(Mat,IS,IS,double); 1008 extern int MatIncreaseOverlap_SeqBAIJ(Mat,int,IS*,int); 1009 extern int MatGetSubMatrix_SeqBAIJ(Mat,IS,IS,int,MatGetSubMatrixCall,Mat*); 1010 extern int MatGetSubMatrices_SeqBAIJ(Mat,int,IS*,IS*,MatGetSubMatrixCall,Mat**); 1011 extern int MatMultTrans_SeqBAIJ(Mat,Vec,Vec); 1012 extern int MatMultTransAdd_SeqBAIJ(Mat,Vec,Vec,Vec); 1013 extern int MatScale_SeqBAIJ(Scalar*,Mat); 1014 extern int MatNorm_SeqBAIJ(Mat,NormType,double *); 1015 extern int MatEqual_SeqBAIJ(Mat,Mat, PetscTruth*); 1016 extern int MatGetDiagonal_SeqBAIJ(Mat,Vec); 1017 extern int MatDiagonalScale_SeqBAIJ(Mat,Vec,Vec); 1018 extern int MatGetInfo_SeqBAIJ(Mat,MatInfoType,MatInfo *); 1019 extern int MatZeroEntries_SeqBAIJ(Mat); 1020 1021 extern int MatSolve_SeqBAIJ_N(Mat,Vec,Vec); 1022 extern int MatSolve_SeqBAIJ_1(Mat,Vec,Vec); 1023 extern int MatSolve_SeqBAIJ_2(Mat,Vec,Vec); 1024 extern int MatSolve_SeqBAIJ_3(Mat,Vec,Vec); 1025 extern int MatSolve_SeqBAIJ_4(Mat,Vec,Vec); 1026 extern int MatSolve_SeqBAIJ_5(Mat,Vec,Vec); 1027 extern int MatSolve_SeqBAIJ_7(Mat,Vec,Vec); 1028 1029 extern int MatLUFactorNumeric_SeqBAIJ_N(Mat,Mat*); 1030 extern int MatLUFactorNumeric_SeqBAIJ_1(Mat,Mat*); 1031 extern int MatLUFactorNumeric_SeqBAIJ_2(Mat,Mat*); 1032 extern int MatLUFactorNumeric_SeqBAIJ_3(Mat,Mat*); 1033 extern int MatLUFactorNumeric_SeqBAIJ_4(Mat,Mat*); 1034 extern int MatLUFactorNumeric_SeqBAIJ_5(Mat,Mat*); 1035 extern int MatSolve_SeqBAIJ_4_NaturalOrdering(Mat,Vec,Vec); 1036 1037 extern int MatMult_SeqBAIJ_1(Mat,Vec,Vec); 1038 extern int MatMult_SeqBAIJ_2(Mat,Vec,Vec); 1039 extern int MatMult_SeqBAIJ_3(Mat,Vec,Vec); 1040 extern int MatMult_SeqBAIJ_4(Mat,Vec,Vec); 1041 extern int MatMult_SeqBAIJ_5(Mat,Vec,Vec); 1042 extern int MatMult_SeqBAIJ_7(Mat,Vec,Vec); 1043 extern int MatMult_SeqBAIJ_N(Mat,Vec,Vec); 1044 1045 extern int MatMultAdd_SeqBAIJ_1(Mat,Vec,Vec,Vec); 1046 extern int MatMultAdd_SeqBAIJ_2(Mat,Vec,Vec,Vec); 1047 extern int MatMultAdd_SeqBAIJ_3(Mat,Vec,Vec,Vec); 1048 extern int MatMultAdd_SeqBAIJ_4(Mat,Vec,Vec,Vec); 1049 extern int MatMultAdd_SeqBAIJ_5(Mat,Vec,Vec,Vec); 1050 extern int MatMultAdd_SeqBAIJ_7(Mat,Vec,Vec,Vec); 1051 extern int MatMultAdd_SeqBAIJ_N(Mat,Vec,Vec,Vec); 1052 1053 /* 1054 note: This can only work for identity for row and col. It would 1055 be good to check this and otherwise generate an error. 1056 */ 1057 #undef __FUNC__ 1058 #define __FUNC__ "MatILUFactor_SeqBAIJ" 1059 int MatILUFactor_SeqBAIJ(Mat inA,IS row,IS col,double efill,int fill) 1060 { 1061 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) inA->data; 1062 Mat outA; 1063 int ierr; 1064 1065 PetscFunctionBegin; 1066 if (fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1067 1068 outA = inA; 1069 inA->factor = FACTOR_LU; 1070 a->row = row; 1071 a->col = col; 1072 1073 if (!a->solve_work) { 1074 a->solve_work = (Scalar *) PetscMalloc((a->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1075 PLogObjectMemory(inA,(a->m+a->bs)*sizeof(Scalar)); 1076 } 1077 1078 if (!a->diag) { 1079 ierr = MatMarkDiag_SeqBAIJ(inA); CHKERRQ(ierr); 1080 } 1081 ierr = MatLUFactorNumeric(inA,&outA); CHKERRQ(ierr); 1082 1083 /* 1084 Blocksize 4 has a special faster solver for ILU(0) factorization 1085 with natural ordering 1086 */ 1087 if (a->bs == 4) { 1088 inA->ops.solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 1089 } 1090 1091 PetscFunctionReturn(0); 1092 } 1093 #undef __FUNC__ 1094 #define __FUNC__ "MatPrintHelp_SeqBAIJ" 1095 int MatPrintHelp_SeqBAIJ(Mat A) 1096 { 1097 static int called = 0; 1098 MPI_Comm comm = A->comm; 1099 1100 PetscFunctionBegin; 1101 if (called) {PetscFunctionReturn(0);} else called = 1; 1102 (*PetscHelpPrintf)(comm," Options for MATSEQBAIJ and MATMPIBAIJ matrix formats (the defaults):\n"); 1103 (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n"); 1104 PetscFunctionReturn(0); 1105 } 1106 1107 /* -------------------------------------------------------------------*/ 1108 static struct _MatOps MatOps = {MatSetValues_SeqBAIJ, 1109 MatGetRow_SeqBAIJ,MatRestoreRow_SeqBAIJ, 1110 MatMult_SeqBAIJ_N,MatMultAdd_SeqBAIJ_N, 1111 MatMultTrans_SeqBAIJ,MatMultTransAdd_SeqBAIJ, 1112 MatSolve_SeqBAIJ_N,0, 1113 0,0, 1114 MatLUFactor_SeqBAIJ,0, 1115 0, 1116 MatTranspose_SeqBAIJ, 1117 MatGetInfo_SeqBAIJ,MatEqual_SeqBAIJ, 1118 MatGetDiagonal_SeqBAIJ,MatDiagonalScale_SeqBAIJ,MatNorm_SeqBAIJ, 1119 0,MatAssemblyEnd_SeqBAIJ, 1120 0, 1121 MatSetOption_SeqBAIJ,MatZeroEntries_SeqBAIJ,MatZeroRows_SeqBAIJ, 1122 MatLUFactorSymbolic_SeqBAIJ,MatLUFactorNumeric_SeqBAIJ_N,0,0, 1123 MatGetSize_SeqBAIJ,MatGetSize_SeqBAIJ,MatGetOwnershipRange_SeqBAIJ, 1124 MatILUFactorSymbolic_SeqBAIJ,0, 1125 0,0, 1126 MatConvertSameType_SeqBAIJ,0,0, 1127 MatILUFactor_SeqBAIJ,0,0, 1128 MatGetSubMatrices_SeqBAIJ,MatIncreaseOverlap_SeqBAIJ, 1129 MatGetValues_SeqBAIJ,0, 1130 MatPrintHelp_SeqBAIJ,MatScale_SeqBAIJ, 1131 0,0,0,MatGetBlockSize_SeqBAIJ, 1132 MatGetRowIJ_SeqBAIJ, 1133 MatRestoreRowIJ_SeqBAIJ, 1134 0,0,0,0,0,0, 1135 MatSetValuesBlocked_SeqBAIJ, 1136 MatGetSubMatrix_SeqBAIJ}; 1137 1138 #undef __FUNC__ 1139 #define __FUNC__ "MatCreateSeqBAIJ" 1140 /*@C 1141 MatCreateSeqBAIJ - Creates a sparse matrix in block AIJ (block 1142 compressed row) format. For good matrix assembly performance the 1143 user should preallocate the matrix storage by setting the parameter nz 1144 (or the array nzz). By setting these parameters accurately, performance 1145 during matrix assembly can be increased by more than a factor of 50. 1146 1147 Input Parameters: 1148 . comm - MPI communicator, set to PETSC_COMM_SELF 1149 . bs - size of block 1150 . m - number of rows 1151 . n - number of columns 1152 . nz - number of block nonzeros per block row (same for all rows) 1153 . nzz - array containing the number of block nonzeros in the various block rows 1154 (possibly different for each block row) or PETSC_NULL 1155 1156 Output Parameter: 1157 . A - the matrix 1158 1159 Options Database Keys: 1160 $ -mat_no_unroll - uses code that does not unroll the loops in the 1161 $ block calculations (much slower) 1162 $ -mat_block_size - size of the blocks to use 1163 1164 Notes: 1165 The block AIJ format is fully compatible with standard Fortran 77 1166 storage. That is, the stored row and column indices can begin at 1167 either one (as in Fortran) or zero. See the users' manual for details. 1168 1169 Specify the preallocated storage with either nz or nnz (not both). 1170 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1171 allocation. For additional details, see the users manual chapter on 1172 matrices. 1173 1174 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues() 1175 @*/ 1176 int MatCreateSeqBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz, Mat *A) 1177 { 1178 Mat B; 1179 Mat_SeqBAIJ *b; 1180 int i,len,ierr,flg,mbs=m/bs,nbs=n/bs,bs2=bs*bs,size; 1181 1182 PetscFunctionBegin; 1183 MPI_Comm_size(comm,&size); 1184 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1185 1186 ierr = OptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,&flg);CHKERRQ(ierr); 1187 1188 if (mbs*bs!=m || nbs*bs!=n) { 1189 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Number rows, cols must be divisible by blocksize"); 1190 } 1191 1192 *A = 0; 1193 PetscHeaderCreate(B,_p_Mat,MAT_COOKIE,MATSEQBAIJ,comm,MatDestroy,MatView); 1194 PLogObjectCreate(B); 1195 B->data = (void *) (b = PetscNew(Mat_SeqBAIJ)); CHKPTRQ(b); 1196 PetscMemzero(b,sizeof(Mat_SeqBAIJ)); 1197 PetscMemcpy(&B->ops,&MatOps,sizeof(struct _MatOps)); 1198 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg); CHKERRQ(ierr); 1199 if (!flg) { 1200 switch (bs) { 1201 case 1: 1202 B->ops.lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 1203 B->ops.solve = MatSolve_SeqBAIJ_1; 1204 B->ops.mult = MatMult_SeqBAIJ_1; 1205 B->ops.multadd = MatMultAdd_SeqBAIJ_1; 1206 break; 1207 case 2: 1208 B->ops.lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2; 1209 B->ops.solve = MatSolve_SeqBAIJ_2; 1210 B->ops.mult = MatMult_SeqBAIJ_2; 1211 B->ops.multadd = MatMultAdd_SeqBAIJ_2; 1212 break; 1213 case 3: 1214 B->ops.lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3; 1215 B->ops.solve = MatSolve_SeqBAIJ_3; 1216 B->ops.mult = MatMult_SeqBAIJ_3; 1217 B->ops.multadd = MatMultAdd_SeqBAIJ_3; 1218 break; 1219 case 4: 1220 B->ops.lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4; 1221 B->ops.solve = MatSolve_SeqBAIJ_4; 1222 B->ops.mult = MatMult_SeqBAIJ_4; 1223 B->ops.multadd = MatMultAdd_SeqBAIJ_4; 1224 break; 1225 case 5: 1226 B->ops.lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5; 1227 B->ops.solve = MatSolve_SeqBAIJ_5; 1228 B->ops.mult = MatMult_SeqBAIJ_5; 1229 B->ops.multadd = MatMultAdd_SeqBAIJ_5; 1230 break; 1231 case 7: 1232 B->ops.mult = MatMult_SeqBAIJ_7; 1233 B->ops.solve = MatSolve_SeqBAIJ_7; 1234 B->ops.multadd = MatMultAdd_SeqBAIJ_7; 1235 break; 1236 } 1237 } 1238 B->destroy = MatDestroy_SeqBAIJ; 1239 B->view = MatView_SeqBAIJ; 1240 B->factor = 0; 1241 B->lupivotthreshold = 1.0; 1242 B->mapping = 0; 1243 b->row = 0; 1244 b->col = 0; 1245 b->reallocs = 0; 1246 1247 b->m = m; B->m = m; B->M = m; 1248 b->n = n; B->n = n; B->N = n; 1249 b->mbs = mbs; 1250 b->nbs = nbs; 1251 b->imax = (int *) PetscMalloc( (mbs+1)*sizeof(int) ); CHKPTRQ(b->imax); 1252 if (nnz == PETSC_NULL) { 1253 if (nz == PETSC_DEFAULT) nz = 5; 1254 else if (nz <= 0) nz = 1; 1255 for ( i=0; i<mbs; i++ ) b->imax[i] = nz; 1256 nz = nz*mbs; 1257 } else { 1258 nz = 0; 1259 for ( i=0; i<mbs; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1260 } 1261 1262 /* allocate the matrix space */ 1263 len = nz*sizeof(int) + nz*bs2*sizeof(Scalar) + (b->m+1)*sizeof(int); 1264 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1265 PetscMemzero(b->a,nz*bs2*sizeof(Scalar)); 1266 b->j = (int *) (b->a + nz*bs2); 1267 PetscMemzero(b->j,nz*sizeof(int)); 1268 b->i = b->j + nz; 1269 b->singlemalloc = 1; 1270 1271 b->i[0] = 0; 1272 for (i=1; i<mbs+1; i++) { 1273 b->i[i] = b->i[i-1] + b->imax[i-1]; 1274 } 1275 1276 /* b->ilen will count nonzeros in each block row so far. */ 1277 b->ilen = (int *) PetscMalloc((mbs+1)*sizeof(int)); 1278 PLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1279 for ( i=0; i<mbs; i++ ) { b->ilen[i] = 0;} 1280 1281 b->bs = bs; 1282 b->bs2 = bs2; 1283 b->mbs = mbs; 1284 b->nz = 0; 1285 b->maxnz = nz*bs2; 1286 b->sorted = 0; 1287 b->roworiented = 1; 1288 b->nonew = 0; 1289 b->diag = 0; 1290 b->solve_work = 0; 1291 b->mult_work = 0; 1292 b->spptr = 0; 1293 B->info.nz_unneeded = (double)b->maxnz; 1294 1295 *A = B; 1296 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1297 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1298 PetscFunctionReturn(0); 1299 } 1300 1301 #undef __FUNC__ 1302 #define __FUNC__ "MatConvertSameType_SeqBAIJ" 1303 int MatConvertSameType_SeqBAIJ(Mat A,Mat *B,int cpvalues) 1304 { 1305 Mat C; 1306 Mat_SeqBAIJ *c,*a = (Mat_SeqBAIJ *) A->data; 1307 int i,len, mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1308 1309 PetscFunctionBegin; 1310 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,0,"Corrupt matrix"); 1311 1312 *B = 0; 1313 PetscHeaderCreate(C,_p_Mat,MAT_COOKIE,MATSEQBAIJ,A->comm,MatDestroy,MatView); 1314 PLogObjectCreate(C); 1315 C->data = (void *) (c = PetscNew(Mat_SeqBAIJ)); CHKPTRQ(c); 1316 PetscMemcpy(&C->ops,&A->ops,sizeof(struct _MatOps)); 1317 C->destroy = MatDestroy_SeqBAIJ; 1318 C->view = MatView_SeqBAIJ; 1319 C->factor = A->factor; 1320 c->row = 0; 1321 c->col = 0; 1322 C->assembled = PETSC_TRUE; 1323 1324 c->m = C->m = a->m; 1325 c->n = C->n = a->n; 1326 C->M = a->m; 1327 C->N = a->n; 1328 1329 c->bs = a->bs; 1330 c->bs2 = a->bs2; 1331 c->mbs = a->mbs; 1332 c->nbs = a->nbs; 1333 1334 c->imax = (int *) PetscMalloc((mbs+1)*sizeof(int)); CHKPTRQ(c->imax); 1335 c->ilen = (int *) PetscMalloc((mbs+1)*sizeof(int)); CHKPTRQ(c->ilen); 1336 for ( i=0; i<mbs; i++ ) { 1337 c->imax[i] = a->imax[i]; 1338 c->ilen[i] = a->ilen[i]; 1339 } 1340 1341 /* allocate the matrix space */ 1342 c->singlemalloc = 1; 1343 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(Scalar) + sizeof(int)); 1344 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1345 c->j = (int *) (c->a + nz*bs2); 1346 c->i = c->j + nz; 1347 PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int)); 1348 if (mbs > 0) { 1349 PetscMemcpy(c->j,a->j,nz*sizeof(int)); 1350 if (cpvalues == COPY_VALUES) { 1351 PetscMemcpy(c->a,a->a,bs2*nz*sizeof(Scalar)); 1352 } 1353 } 1354 1355 PLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1356 c->sorted = a->sorted; 1357 c->roworiented = a->roworiented; 1358 c->nonew = a->nonew; 1359 1360 if (a->diag) { 1361 c->diag = (int *) PetscMalloc( (mbs+1)*sizeof(int) ); CHKPTRQ(c->diag); 1362 PLogObjectMemory(C,(mbs+1)*sizeof(int)); 1363 for ( i=0; i<mbs; i++ ) { 1364 c->diag[i] = a->diag[i]; 1365 } 1366 } 1367 else c->diag = 0; 1368 c->nz = a->nz; 1369 c->maxnz = a->maxnz; 1370 c->solve_work = 0; 1371 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1372 c->mult_work = 0; 1373 *B = C; 1374 PetscFunctionReturn(0); 1375 } 1376 1377 #undef __FUNC__ 1378 #define __FUNC__ "MatLoad_SeqBAIJ" 1379 int MatLoad_SeqBAIJ(Viewer viewer,MatType type,Mat *A) 1380 { 1381 Mat_SeqBAIJ *a; 1382 Mat B; 1383 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1,flg; 1384 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount; 1385 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1386 int *masked, nmask,tmp,bs2,ishift; 1387 Scalar *aa; 1388 MPI_Comm comm = ((PetscObject) viewer)->comm; 1389 1390 PetscFunctionBegin; 1391 ierr = OptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,&flg);CHKERRQ(ierr); 1392 bs2 = bs*bs; 1393 1394 MPI_Comm_size(comm,&size); 1395 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 1396 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1397 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 1398 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not Mat object"); 1399 M = header[1]; N = header[2]; nz = header[3]; 1400 1401 if (header[3] < 0) { 1402 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format, cannot load as SeqBAIJ"); 1403 } 1404 1405 if (M != N) SETERRQ(PETSC_ERR_SUP,0,"Can only do square matrices"); 1406 1407 /* 1408 This code adds extra rows to make sure the number of rows is 1409 divisible by the blocksize 1410 */ 1411 mbs = M/bs; 1412 extra_rows = bs - M + bs*(mbs); 1413 if (extra_rows == bs) extra_rows = 0; 1414 else mbs++; 1415 if (extra_rows) { 1416 PLogInfo(0,"MatLoad_SeqBAIJ:Padding loaded matrix to match blocksize\n"); 1417 } 1418 1419 /* read in row lengths */ 1420 rowlengths = (int*) PetscMalloc((M+extra_rows)*sizeof(int));CHKPTRQ(rowlengths); 1421 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 1422 for ( i=0; i<extra_rows; i++ ) rowlengths[M+i] = 1; 1423 1424 /* read in column indices */ 1425 jj = (int*) PetscMalloc( (nz+extra_rows)*sizeof(int) ); CHKPTRQ(jj); 1426 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT); CHKERRQ(ierr); 1427 for ( i=0; i<extra_rows; i++ ) jj[nz+i] = M+i; 1428 1429 /* loop over row lengths determining block row lengths */ 1430 browlengths = (int *) PetscMalloc(mbs*sizeof(int));CHKPTRQ(browlengths); 1431 PetscMemzero(browlengths,mbs*sizeof(int)); 1432 mask = (int *) PetscMalloc( 2*mbs*sizeof(int) ); CHKPTRQ(mask); 1433 PetscMemzero(mask,mbs*sizeof(int)); 1434 masked = mask + mbs; 1435 rowcount = 0; nzcount = 0; 1436 for ( i=0; i<mbs; i++ ) { 1437 nmask = 0; 1438 for ( j=0; j<bs; j++ ) { 1439 kmax = rowlengths[rowcount]; 1440 for ( k=0; k<kmax; k++ ) { 1441 tmp = jj[nzcount++]/bs; 1442 if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;} 1443 } 1444 rowcount++; 1445 } 1446 browlengths[i] += nmask; 1447 /* zero out the mask elements we set */ 1448 for ( j=0; j<nmask; j++ ) mask[masked[j]] = 0; 1449 } 1450 1451 /* create our matrix */ 1452 ierr = MatCreateSeqBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1453 B = *A; 1454 a = (Mat_SeqBAIJ *) B->data; 1455 1456 /* set matrix "i" values */ 1457 a->i[0] = 0; 1458 for ( i=1; i<= mbs; i++ ) { 1459 a->i[i] = a->i[i-1] + browlengths[i-1]; 1460 a->ilen[i-1] = browlengths[i-1]; 1461 } 1462 a->nz = 0; 1463 for ( i=0; i<mbs; i++ ) a->nz += browlengths[i]; 1464 1465 /* read in nonzero values */ 1466 aa = (Scalar *) PetscMalloc((nz+extra_rows)*sizeof(Scalar));CHKPTRQ(aa); 1467 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR); CHKERRQ(ierr); 1468 for ( i=0; i<extra_rows; i++ ) aa[nz+i] = 1.0; 1469 1470 /* set "a" and "j" values into matrix */ 1471 nzcount = 0; jcount = 0; 1472 for ( i=0; i<mbs; i++ ) { 1473 nzcountb = nzcount; 1474 nmask = 0; 1475 for ( j=0; j<bs; j++ ) { 1476 kmax = rowlengths[i*bs+j]; 1477 for ( k=0; k<kmax; k++ ) { 1478 tmp = jj[nzcount++]/bs; 1479 if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;} 1480 } 1481 rowcount++; 1482 } 1483 /* sort the masked values */ 1484 PetscSortInt(nmask,masked); 1485 1486 /* set "j" values into matrix */ 1487 maskcount = 1; 1488 for ( j=0; j<nmask; j++ ) { 1489 a->j[jcount++] = masked[j]; 1490 mask[masked[j]] = maskcount++; 1491 } 1492 /* set "a" values into matrix */ 1493 ishift = bs2*a->i[i]; 1494 for ( j=0; j<bs; j++ ) { 1495 kmax = rowlengths[i*bs+j]; 1496 for ( k=0; k<kmax; k++ ) { 1497 tmp = jj[nzcountb]/bs ; 1498 block = mask[tmp] - 1; 1499 point = jj[nzcountb] - bs*tmp; 1500 idx = ishift + bs2*block + j + bs*point; 1501 a->a[idx] = aa[nzcountb++]; 1502 } 1503 } 1504 /* zero out the mask elements we set */ 1505 for ( j=0; j<nmask; j++ ) mask[masked[j]] = 0; 1506 } 1507 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Bad binary matrix"); 1508 1509 PetscFree(rowlengths); 1510 PetscFree(browlengths); 1511 PetscFree(aa); 1512 PetscFree(jj); 1513 PetscFree(mask); 1514 1515 B->assembled = PETSC_TRUE; 1516 1517 ierr = MatView_Private(B); CHKERRQ(ierr); 1518 PetscFunctionReturn(0); 1519 } 1520 1521 1522 1523