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