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