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