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