1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: baij.c,v 1.146 1998/10/19 22:18:22 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 1073 extern int MatMult_SeqBAIJ_1(Mat,Vec,Vec); 1074 extern int MatMult_SeqBAIJ_2(Mat,Vec,Vec); 1075 extern int MatMult_SeqBAIJ_3(Mat,Vec,Vec); 1076 extern int MatMult_SeqBAIJ_4(Mat,Vec,Vec); 1077 extern int MatMult_SeqBAIJ_5(Mat,Vec,Vec); 1078 extern int MatMult_SeqBAIJ_7(Mat,Vec,Vec); 1079 extern int MatMult_SeqBAIJ_N(Mat,Vec,Vec); 1080 1081 extern int MatMultAdd_SeqBAIJ_1(Mat,Vec,Vec,Vec); 1082 extern int MatMultAdd_SeqBAIJ_2(Mat,Vec,Vec,Vec); 1083 extern int MatMultAdd_SeqBAIJ_3(Mat,Vec,Vec,Vec); 1084 extern int MatMultAdd_SeqBAIJ_4(Mat,Vec,Vec,Vec); 1085 extern int MatMultAdd_SeqBAIJ_5(Mat,Vec,Vec,Vec); 1086 extern int MatMultAdd_SeqBAIJ_7(Mat,Vec,Vec,Vec); 1087 extern int MatMultAdd_SeqBAIJ_N(Mat,Vec,Vec,Vec); 1088 1089 #undef __FUNC__ 1090 #define __FUNC__ "MatILUFactor_SeqBAIJ" 1091 int MatILUFactor_SeqBAIJ(Mat inA,IS row,IS col,double efill,int fill) 1092 { 1093 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) inA->data; 1094 Mat outA; 1095 int ierr; 1096 PetscTruth row_identity, col_identity; 1097 1098 PetscFunctionBegin; 1099 if (fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1100 ierr = ISIdentity(row,&row_identity); CHKERRQ(ierr); 1101 ierr = ISIdentity(col,&col_identity); CHKERRQ(ierr); 1102 if (!row_identity || !col_identity) { 1103 SETERRQ(1,1,"Row and column permutations must be identity"); 1104 } 1105 1106 outA = inA; 1107 inA->factor = FACTOR_LU; 1108 a->row = row; 1109 a->col = col; 1110 1111 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1112 ierr = ISInvertPermutation(col,&(a->icol)); CHKERRQ(ierr); 1113 PLogObjectParent(inA,a->icol); 1114 1115 if (!a->solve_work) { 1116 a->solve_work = (Scalar *) PetscMalloc((a->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1117 PLogObjectMemory(inA,(a->m+a->bs)*sizeof(Scalar)); 1118 } 1119 1120 if (!a->diag) { 1121 ierr = MatMarkDiag_SeqBAIJ(inA); CHKERRQ(ierr); 1122 } 1123 /* 1124 Blocksize 4 and 5 have a special faster factorization/solver for ILU(0) factorization 1125 with natural ordering 1126 */ 1127 if (a->bs == 4) { 1128 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 1129 inA->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 1130 PLogInfo(inA,"Using special in-place natural ordering factor and solve BS=4\n"); 1131 } else if (a->bs == 5) { 1132 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering; 1133 inA->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering; 1134 PLogInfo(inA,"Using special in-place natural ordering factor and solve BS=5\n"); 1135 } 1136 1137 ierr = MatLUFactorNumeric(inA,&outA); CHKERRQ(ierr); 1138 1139 1140 PetscFunctionReturn(0); 1141 } 1142 #undef __FUNC__ 1143 #define __FUNC__ "MatPrintHelp_SeqBAIJ" 1144 int MatPrintHelp_SeqBAIJ(Mat A) 1145 { 1146 static int called = 0; 1147 MPI_Comm comm = A->comm; 1148 1149 PetscFunctionBegin; 1150 if (called) {PetscFunctionReturn(0);} else called = 1; 1151 (*PetscHelpPrintf)(comm," Options for MATSEQBAIJ and MATMPIBAIJ matrix formats (the defaults):\n"); 1152 (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n"); 1153 PetscFunctionReturn(0); 1154 } 1155 1156 EXTERN_C_BEGIN 1157 #undef __FUNC__ 1158 #define __FUNC__ "MatSeqBAIJSetColumnIndices_SeqBAIJ" 1159 int MatSeqBAIJSetColumnIndices_SeqBAIJ(Mat mat,int *indices) 1160 { 1161 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data; 1162 int i,nz,n; 1163 1164 PetscFunctionBegin; 1165 nz = baij->maxnz; 1166 n = baij->n; 1167 for (i=0; i<nz; i++) { 1168 baij->j[i] = indices[i]; 1169 } 1170 baij->nz = nz; 1171 for ( i=0; i<n; i++ ) { 1172 baij->ilen[i] = baij->imax[i]; 1173 } 1174 1175 PetscFunctionReturn(0); 1176 } 1177 EXTERN_C_END 1178 1179 #undef __FUNC__ 1180 #define __FUNC__ "MatSeqBAIJSetColumnIndices" 1181 /*@ 1182 MatSeqBAIJSetColumnIndices - Set the column indices for all the rows 1183 in the matrix. 1184 1185 Input Parameters: 1186 + mat - the SeqBAIJ matrix 1187 - indices - the column indices 1188 1189 Notes: 1190 This can be called if you have precomputed the nonzero structure of the 1191 matrix and want to provide it to the matrix object to improve the performance 1192 of the MatSetValues() operation. 1193 1194 You MUST have set the correct numbers of nonzeros per row in the call to 1195 MatCreateSeqBAIJ(). 1196 1197 MUST be called before any calls to MatSetValues(); 1198 1199 @*/ 1200 int MatSeqBAIJSetColumnIndices(Mat mat,int *indices) 1201 { 1202 int ierr,(*f)(Mat,int *); 1203 1204 PetscFunctionBegin; 1205 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1206 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqBAIJSetColumnIndices_C",(void **)&f); 1207 CHKERRQ(ierr); 1208 if (f) { 1209 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1210 } else { 1211 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1212 } 1213 PetscFunctionReturn(0); 1214 } 1215 1216 /* -------------------------------------------------------------------*/ 1217 static struct _MatOps MatOps_Values = {MatSetValues_SeqBAIJ, 1218 MatGetRow_SeqBAIJ, 1219 MatRestoreRow_SeqBAIJ, 1220 MatMult_SeqBAIJ_N, 1221 MatMultAdd_SeqBAIJ_N, 1222 MatMultTrans_SeqBAIJ, 1223 MatMultTransAdd_SeqBAIJ, 1224 MatSolve_SeqBAIJ_N, 1225 0, 1226 0, 1227 0, 1228 MatLUFactor_SeqBAIJ, 1229 0, 1230 0, 1231 MatTranspose_SeqBAIJ, 1232 MatGetInfo_SeqBAIJ, 1233 MatEqual_SeqBAIJ, 1234 MatGetDiagonal_SeqBAIJ, 1235 MatDiagonalScale_SeqBAIJ, 1236 MatNorm_SeqBAIJ, 1237 0, 1238 MatAssemblyEnd_SeqBAIJ, 1239 0, 1240 MatSetOption_SeqBAIJ, 1241 MatZeroEntries_SeqBAIJ, 1242 MatZeroRows_SeqBAIJ, 1243 MatLUFactorSymbolic_SeqBAIJ, 1244 MatLUFactorNumeric_SeqBAIJ_N, 1245 0, 1246 0, 1247 MatGetSize_SeqBAIJ, 1248 MatGetSize_SeqBAIJ, 1249 MatGetOwnershipRange_SeqBAIJ, 1250 MatILUFactorSymbolic_SeqBAIJ, 1251 0, 1252 0, 1253 0, 1254 MatDuplicate_SeqBAIJ, 1255 0, 1256 0, 1257 MatILUFactor_SeqBAIJ, 1258 0, 1259 0, 1260 MatGetSubMatrices_SeqBAIJ, 1261 MatIncreaseOverlap_SeqBAIJ, 1262 MatGetValues_SeqBAIJ, 1263 0, 1264 MatPrintHelp_SeqBAIJ, 1265 MatScale_SeqBAIJ, 1266 0, 1267 0, 1268 0, 1269 MatGetBlockSize_SeqBAIJ, 1270 MatGetRowIJ_SeqBAIJ, 1271 MatRestoreRowIJ_SeqBAIJ, 1272 0, 1273 0, 1274 0, 1275 0, 1276 0, 1277 0, 1278 MatSetValuesBlocked_SeqBAIJ, 1279 MatGetSubMatrix_SeqBAIJ, 1280 0, 1281 0, 1282 MatGetMaps_Petsc}; 1283 1284 #undef __FUNC__ 1285 #define __FUNC__ "MatCreateSeqBAIJ" 1286 /*@C 1287 MatCreateSeqBAIJ - Creates a sparse matrix in block AIJ (block 1288 compressed row) format. For good matrix assembly performance the 1289 user should preallocate the matrix storage by setting the parameter nz 1290 (or the array nzz). By setting these parameters accurately, performance 1291 during matrix assembly can be increased by more than a factor of 50. 1292 1293 Collective on MPI_Comm 1294 1295 Input Parameters: 1296 + comm - MPI communicator, set to PETSC_COMM_SELF 1297 . bs - size of block 1298 . m - number of rows 1299 . n - number of columns 1300 . nz - number of block nonzeros per block row (same for all rows) 1301 - nzz - array containing the number of block nonzeros in the various block rows 1302 (possibly different for each block row) or PETSC_NULL 1303 1304 Output Parameter: 1305 . A - the matrix 1306 1307 Options Database Keys: 1308 . -mat_no_unroll - uses code that does not unroll the loops in the 1309 block calculations (much slower) 1310 . -mat_block_size - size of the blocks to use 1311 1312 Notes: 1313 The block AIJ format is fully compatible with standard Fortran 77 1314 storage. That is, the stored row and column indices can begin at 1315 either one (as in Fortran) or zero. See the users' manual for details. 1316 1317 Specify the preallocated storage with either nz or nnz (not both). 1318 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1319 allocation. For additional details, see the users manual chapter on 1320 matrices. 1321 1322 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1323 @*/ 1324 int MatCreateSeqBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz, Mat *A) 1325 { 1326 Mat B; 1327 Mat_SeqBAIJ *b; 1328 int i,len,ierr,flg,mbs=m/bs,nbs=n/bs,bs2=bs*bs,size; 1329 1330 PetscFunctionBegin; 1331 MPI_Comm_size(comm,&size); 1332 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1333 1334 ierr = OptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,&flg);CHKERRQ(ierr); 1335 1336 if (mbs*bs!=m || nbs*bs!=n) { 1337 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Number rows, cols must be divisible by blocksize"); 1338 } 1339 1340 *A = 0; 1341 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQBAIJ,comm,MatDestroy,MatView); 1342 PLogObjectCreate(B); 1343 B->data = (void *) (b = PetscNew(Mat_SeqBAIJ)); CHKPTRQ(b); 1344 PetscMemzero(b,sizeof(Mat_SeqBAIJ)); 1345 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 1346 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg); CHKERRQ(ierr); 1347 if (!flg) { 1348 switch (bs) { 1349 case 1: 1350 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 1351 B->ops->solve = MatSolve_SeqBAIJ_1; 1352 B->ops->mult = MatMult_SeqBAIJ_1; 1353 B->ops->multadd = MatMultAdd_SeqBAIJ_1; 1354 break; 1355 case 2: 1356 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2; 1357 B->ops->solve = MatSolve_SeqBAIJ_2; 1358 B->ops->mult = MatMult_SeqBAIJ_2; 1359 B->ops->multadd = MatMultAdd_SeqBAIJ_2; 1360 break; 1361 case 3: 1362 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3; 1363 B->ops->solve = MatSolve_SeqBAIJ_3; 1364 B->ops->mult = MatMult_SeqBAIJ_3; 1365 B->ops->multadd = MatMultAdd_SeqBAIJ_3; 1366 break; 1367 case 4: 1368 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4; 1369 B->ops->solve = MatSolve_SeqBAIJ_4; 1370 B->ops->mult = MatMult_SeqBAIJ_4; 1371 B->ops->multadd = MatMultAdd_SeqBAIJ_4; 1372 break; 1373 case 5: 1374 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5; 1375 B->ops->solve = MatSolve_SeqBAIJ_5; 1376 B->ops->mult = MatMult_SeqBAIJ_5; 1377 B->ops->multadd = MatMultAdd_SeqBAIJ_5; 1378 break; 1379 case 7: 1380 B->ops->mult = MatMult_SeqBAIJ_7; 1381 B->ops->solve = MatSolve_SeqBAIJ_7; 1382 B->ops->multadd = MatMultAdd_SeqBAIJ_7; 1383 break; 1384 } 1385 } 1386 B->ops->destroy = MatDestroy_SeqBAIJ; 1387 B->ops->view = MatView_SeqBAIJ; 1388 B->factor = 0; 1389 B->lupivotthreshold = 1.0; 1390 B->mapping = 0; 1391 b->row = 0; 1392 b->col = 0; 1393 b->icol = 0; 1394 b->reallocs = 0; 1395 1396 b->m = m; B->m = m; B->M = m; 1397 b->n = n; B->n = n; B->N = n; 1398 1399 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1400 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 1401 1402 b->mbs = mbs; 1403 b->nbs = nbs; 1404 b->imax = (int *) PetscMalloc( (mbs+1)*sizeof(int) ); CHKPTRQ(b->imax); 1405 if (nnz == PETSC_NULL) { 1406 if (nz == PETSC_DEFAULT) nz = 5; 1407 else if (nz <= 0) nz = 1; 1408 for ( i=0; i<mbs; i++ ) b->imax[i] = nz; 1409 nz = nz*mbs; 1410 } else { 1411 nz = 0; 1412 for ( i=0; i<mbs; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1413 } 1414 1415 /* allocate the matrix space */ 1416 len = nz*sizeof(int) + nz*bs2*sizeof(Scalar) + (b->m+1)*sizeof(int); 1417 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1418 PetscMemzero(b->a,nz*bs2*sizeof(Scalar)); 1419 b->j = (int *) (b->a + nz*bs2); 1420 PetscMemzero(b->j,nz*sizeof(int)); 1421 b->i = b->j + nz; 1422 b->singlemalloc = 1; 1423 1424 b->i[0] = 0; 1425 for (i=1; i<mbs+1; i++) { 1426 b->i[i] = b->i[i-1] + b->imax[i-1]; 1427 } 1428 1429 /* b->ilen will count nonzeros in each block row so far. */ 1430 b->ilen = (int *) PetscMalloc((mbs+1)*sizeof(int)); 1431 PLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1432 for ( i=0; i<mbs; i++ ) { b->ilen[i] = 0;} 1433 1434 b->bs = bs; 1435 b->bs2 = bs2; 1436 b->mbs = mbs; 1437 b->nz = 0; 1438 b->maxnz = nz*bs2; 1439 b->sorted = 0; 1440 b->roworiented = 1; 1441 b->nonew = 0; 1442 b->diag = 0; 1443 b->solve_work = 0; 1444 b->mult_work = 0; 1445 b->spptr = 0; 1446 B->info.nz_unneeded = (double)b->maxnz; 1447 1448 *A = B; 1449 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1450 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1451 1452 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqBAIJSetColumnIndices_C", 1453 "MatSeqBAIJSetColumnIndices_SeqBAIJ", 1454 (void*)MatSeqBAIJSetColumnIndices_SeqBAIJ);CHKERRQ(ierr); 1455 1456 PetscFunctionReturn(0); 1457 } 1458 1459 #undef __FUNC__ 1460 #define __FUNC__ "MatDuplicate_SeqBAIJ" 1461 int MatDuplicate_SeqBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1462 { 1463 Mat C; 1464 Mat_SeqBAIJ *c,*a = (Mat_SeqBAIJ *) A->data; 1465 int i,len, mbs = a->mbs,nz = a->nz,bs2 =a->bs2,ierr; 1466 1467 PetscFunctionBegin; 1468 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,0,"Corrupt matrix"); 1469 1470 *B = 0; 1471 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQBAIJ,A->comm,MatDestroy,MatView); 1472 PLogObjectCreate(C); 1473 C->data = (void *) (c = PetscNew(Mat_SeqBAIJ)); CHKPTRQ(c); 1474 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 1475 C->ops->destroy = MatDestroy_SeqBAIJ; 1476 C->ops->view = MatView_SeqBAIJ; 1477 C->factor = A->factor; 1478 c->row = 0; 1479 c->col = 0; 1480 C->assembled = PETSC_TRUE; 1481 1482 c->m = C->m = a->m; 1483 c->n = C->n = a->n; 1484 C->M = a->m; 1485 C->N = a->n; 1486 1487 c->bs = a->bs; 1488 c->bs2 = a->bs2; 1489 c->mbs = a->mbs; 1490 c->nbs = a->nbs; 1491 1492 c->imax = (int *) PetscMalloc((mbs+1)*sizeof(int)); CHKPTRQ(c->imax); 1493 c->ilen = (int *) PetscMalloc((mbs+1)*sizeof(int)); CHKPTRQ(c->ilen); 1494 for ( i=0; i<mbs; i++ ) { 1495 c->imax[i] = a->imax[i]; 1496 c->ilen[i] = a->ilen[i]; 1497 } 1498 1499 /* allocate the matrix space */ 1500 c->singlemalloc = 1; 1501 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(Scalar) + sizeof(int)); 1502 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1503 c->j = (int *) (c->a + nz*bs2); 1504 c->i = c->j + nz; 1505 PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int)); 1506 if (mbs > 0) { 1507 PetscMemcpy(c->j,a->j,nz*sizeof(int)); 1508 if (cpvalues == MAT_COPY_VALUES) { 1509 PetscMemcpy(c->a,a->a,bs2*nz*sizeof(Scalar)); 1510 } else { 1511 PetscMemzero(c->a,bs2*nz*sizeof(Scalar)); 1512 } 1513 } 1514 1515 PLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1516 c->sorted = a->sorted; 1517 c->roworiented = a->roworiented; 1518 c->nonew = a->nonew; 1519 1520 if (a->diag) { 1521 c->diag = (int *) PetscMalloc( (mbs+1)*sizeof(int) ); CHKPTRQ(c->diag); 1522 PLogObjectMemory(C,(mbs+1)*sizeof(int)); 1523 for ( i=0; i<mbs; i++ ) { 1524 c->diag[i] = a->diag[i]; 1525 } 1526 } 1527 else c->diag = 0; 1528 c->nz = a->nz; 1529 c->maxnz = a->maxnz; 1530 c->solve_work = 0; 1531 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1532 c->mult_work = 0; 1533 *B = C; 1534 ierr = PetscObjectComposeFunction((PetscObject)C,"MatSeqBAIJSetColumnIndices_C", 1535 "MatSeqBAIJSetColumnIndices_SeqBAIJ", 1536 (void*)MatSeqBAIJSetColumnIndices_SeqBAIJ);CHKERRQ(ierr); 1537 PetscFunctionReturn(0); 1538 } 1539 1540 #undef __FUNC__ 1541 #define __FUNC__ "MatLoad_SeqBAIJ" 1542 int MatLoad_SeqBAIJ(Viewer viewer,MatType type,Mat *A) 1543 { 1544 Mat_SeqBAIJ *a; 1545 Mat B; 1546 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1,flg; 1547 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount; 1548 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1549 int *masked, nmask,tmp,bs2,ishift; 1550 Scalar *aa; 1551 MPI_Comm comm = ((PetscObject) viewer)->comm; 1552 1553 PetscFunctionBegin; 1554 ierr = OptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,&flg);CHKERRQ(ierr); 1555 bs2 = bs*bs; 1556 1557 MPI_Comm_size(comm,&size); 1558 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 1559 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1560 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 1561 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not Mat object"); 1562 M = header[1]; N = header[2]; nz = header[3]; 1563 1564 if (header[3] < 0) { 1565 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format, cannot load as SeqBAIJ"); 1566 } 1567 1568 if (M != N) SETERRQ(PETSC_ERR_SUP,0,"Can only do square matrices"); 1569 1570 /* 1571 This code adds extra rows to make sure the number of rows is 1572 divisible by the blocksize 1573 */ 1574 mbs = M/bs; 1575 extra_rows = bs - M + bs*(mbs); 1576 if (extra_rows == bs) extra_rows = 0; 1577 else mbs++; 1578 if (extra_rows) { 1579 PLogInfo(0,"MatLoad_SeqBAIJ:Padding loaded matrix to match blocksize\n"); 1580 } 1581 1582 /* read in row lengths */ 1583 rowlengths = (int*) PetscMalloc((M+extra_rows)*sizeof(int));CHKPTRQ(rowlengths); 1584 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 1585 for ( i=0; i<extra_rows; i++ ) rowlengths[M+i] = 1; 1586 1587 /* read in column indices */ 1588 jj = (int*) PetscMalloc( (nz+extra_rows)*sizeof(int) ); CHKPTRQ(jj); 1589 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT); CHKERRQ(ierr); 1590 for ( i=0; i<extra_rows; i++ ) jj[nz+i] = M+i; 1591 1592 /* loop over row lengths determining block row lengths */ 1593 browlengths = (int *) PetscMalloc(mbs*sizeof(int));CHKPTRQ(browlengths); 1594 PetscMemzero(browlengths,mbs*sizeof(int)); 1595 mask = (int *) PetscMalloc( 2*mbs*sizeof(int) ); CHKPTRQ(mask); 1596 PetscMemzero(mask,mbs*sizeof(int)); 1597 masked = mask + mbs; 1598 rowcount = 0; nzcount = 0; 1599 for ( i=0; i<mbs; i++ ) { 1600 nmask = 0; 1601 for ( j=0; j<bs; j++ ) { 1602 kmax = rowlengths[rowcount]; 1603 for ( k=0; k<kmax; k++ ) { 1604 tmp = jj[nzcount++]/bs; 1605 if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;} 1606 } 1607 rowcount++; 1608 } 1609 browlengths[i] += nmask; 1610 /* zero out the mask elements we set */ 1611 for ( j=0; j<nmask; j++ ) mask[masked[j]] = 0; 1612 } 1613 1614 /* create our matrix */ 1615 ierr = MatCreateSeqBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1616 B = *A; 1617 a = (Mat_SeqBAIJ *) B->data; 1618 1619 /* set matrix "i" values */ 1620 a->i[0] = 0; 1621 for ( i=1; i<= mbs; i++ ) { 1622 a->i[i] = a->i[i-1] + browlengths[i-1]; 1623 a->ilen[i-1] = browlengths[i-1]; 1624 } 1625 a->nz = 0; 1626 for ( i=0; i<mbs; i++ ) a->nz += browlengths[i]; 1627 1628 /* read in nonzero values */ 1629 aa = (Scalar *) PetscMalloc((nz+extra_rows)*sizeof(Scalar));CHKPTRQ(aa); 1630 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR); CHKERRQ(ierr); 1631 for ( i=0; i<extra_rows; i++ ) aa[nz+i] = 1.0; 1632 1633 /* set "a" and "j" values into matrix */ 1634 nzcount = 0; jcount = 0; 1635 for ( i=0; i<mbs; i++ ) { 1636 nzcountb = nzcount; 1637 nmask = 0; 1638 for ( j=0; j<bs; j++ ) { 1639 kmax = rowlengths[i*bs+j]; 1640 for ( k=0; k<kmax; k++ ) { 1641 tmp = jj[nzcount++]/bs; 1642 if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;} 1643 } 1644 rowcount++; 1645 } 1646 /* sort the masked values */ 1647 PetscSortInt(nmask,masked); 1648 1649 /* set "j" values into matrix */ 1650 maskcount = 1; 1651 for ( j=0; j<nmask; j++ ) { 1652 a->j[jcount++] = masked[j]; 1653 mask[masked[j]] = maskcount++; 1654 } 1655 /* set "a" values into matrix */ 1656 ishift = bs2*a->i[i]; 1657 for ( j=0; j<bs; j++ ) { 1658 kmax = rowlengths[i*bs+j]; 1659 for ( k=0; k<kmax; k++ ) { 1660 tmp = jj[nzcountb]/bs ; 1661 block = mask[tmp] - 1; 1662 point = jj[nzcountb] - bs*tmp; 1663 idx = ishift + bs2*block + j + bs*point; 1664 a->a[idx] = aa[nzcountb++]; 1665 } 1666 } 1667 /* zero out the mask elements we set */ 1668 for ( j=0; j<nmask; j++ ) mask[masked[j]] = 0; 1669 } 1670 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Bad binary matrix"); 1671 1672 PetscFree(rowlengths); 1673 PetscFree(browlengths); 1674 PetscFree(aa); 1675 PetscFree(jj); 1676 PetscFree(mask); 1677 1678 B->assembled = PETSC_TRUE; 1679 1680 ierr = MatView_Private(B); CHKERRQ(ierr); 1681 PetscFunctionReturn(0); 1682 } 1683 1684 1685 1686