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