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