1 /*$Id: sbaij.c,v 1.44 2000/11/01 20:41:07 hzhang 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; 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_SYMMETRIC || 153 op == MAT_STRUCTURALLY_SYMMETRIC || 154 op == MAT_YES_NEW_DIAGONALS || 155 op == MAT_IGNORE_OFF_PROC_ENTRIES || 156 op == MAT_USE_HASH_TABLE) { 157 PetscLogInfo(A,"MatSetOption_SeqSBAIJ:Option ignored\n"); 158 } else if (op == MAT_NO_NEW_DIAGONALS) { 159 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 160 } else { 161 SETERRQ(PETSC_ERR_SUP,"unknown option"); 162 } 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNC__ 167 #define __FUNC__ "MatGetOwnershipRange_SeqSBAIJ" 168 int MatGetOwnershipRange_SeqSBAIJ(Mat A,int *m,int *n) 169 { 170 PetscFunctionBegin; 171 *m = 0; 172 *n = A->m; 173 PetscFunctionReturn(0); 174 } 175 176 #undef __FUNC__ 177 #define __FUNC__ "MatGetRow_SeqSBAIJ" 178 int MatGetRow_SeqSBAIJ(Mat A,int row,int *ncols,int **cols,Scalar **v) 179 { 180 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 181 int itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*cols_i,bs2; 182 MatScalar *aa,*aa_i; 183 Scalar *v_i; 184 185 PetscFunctionBegin; 186 bs = a->bs; 187 ai = a->i; 188 aj = a->j; 189 aa = a->a; 190 bs2 = a->bs2; 191 192 if (row < 0 || row >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row out of range"); 193 194 bn = row/bs; /* Block number */ 195 bp = row % bs; /* Block position */ 196 M = ai[bn+1] - ai[bn]; 197 *ncols = bs*M; 198 199 if (v) { 200 *v = 0; 201 if (*ncols) { 202 ierr = PetscMalloc((*ncols+row)*sizeof(Scalar),& *v );CHKERRQ(ierr); 203 for (i=0; i<M; i++) { /* for each block in the block row */ 204 v_i = *v + i*bs; 205 aa_i = aa + bs2*(ai[bn] + i); 206 for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];} 207 } 208 } 209 } 210 211 if (cols) { 212 *cols = 0; 213 if (*ncols) { 214 ierr = PetscMalloc((*ncols+row)*sizeof(int),& *cols );CHKERRQ(ierr); 215 for (i=0; i<M; i++) { /* for each block in the block row */ 216 cols_i = *cols + i*bs; 217 itmp = bs*aj[ai[bn] + i]; 218 for (j=0; j<bs; j++) {cols_i[j] = itmp++;} 219 } 220 } 221 } 222 223 /*search column A(0:row-1,row) (=A(row,0:row-1)). Could be expensive! */ 224 /* this segment is currently removed, so only entries in the upper triangle are obtained */ 225 #ifdef column_search 226 v_i = *v + M*bs; 227 cols_i = *cols + M*bs; 228 for (i=0; i<bn; i++){ /* for each block row */ 229 M = ai[i+1] - ai[i]; 230 for (j=0; j<M; j++){ 231 itmp = aj[ai[i] + j]; /* block column value */ 232 if (itmp == bn){ 233 aa_i = aa + bs2*(ai[i] + j) + bs*bp; 234 for (k=0; k<bs; k++) { 235 *cols_i++ = i*bs+k; 236 *v_i++ = aa_i[k]; 237 } 238 *ncols += bs; 239 break; 240 } 241 } 242 } 243 #endif 244 245 PetscFunctionReturn(0); 246 } 247 248 #undef __FUNC__ 249 #define __FUNC__ "MatRestoreRow_SeqSBAIJ" 250 int MatRestoreRow_SeqSBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 251 { 252 int ierr; 253 254 PetscFunctionBegin; 255 if (idx) {if (*idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 256 if (v) {if (*v) {ierr = PetscFree(*v);CHKERRQ(ierr);}} 257 PetscFunctionReturn(0); 258 } 259 260 #undef __FUNC__ 261 #define __FUNC__ "MatTranspose_SeqSBAIJ" 262 int MatTranspose_SeqSBAIJ(Mat A,Mat *B) 263 { 264 PetscFunctionBegin; 265 SETERRQ(1,"Matrix is symmetric. MatTranspose() should not be called"); 266 /* PetscFunctionReturn(0); */ 267 } 268 269 #undef __FUNC__ 270 #define __FUNC__ "MatView_SeqSBAIJ_Binary" 271 static int MatView_SeqSBAIJ_Binary(Mat A,PetscViewer viewer) 272 { 273 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 274 int i,fd,*col_lens,ierr,bs = a->bs,count,*jj,j,k,l,bs2=a->bs2; 275 Scalar *aa; 276 FILE *file; 277 278 PetscFunctionBegin; 279 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 280 col_lens = (int*)PetscMalloc((4+A->m)*sizeof(int));CHKERRQ(ierr); 281 col_lens[0] = MAT_COOKIE; 282 283 col_lens[1] = A->m; 284 col_lens[2] = A->m; 285 col_lens[3] = a->s_nz*bs2; 286 287 /* store lengths of each row and write (including header) to file */ 288 count = 0; 289 for (i=0; i<a->mbs; i++) { 290 for (j=0; j<bs; j++) { 291 col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]); 292 } 293 } 294 ierr = PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,1);CHKERRQ(ierr); 295 ierr = PetscFree(col_lens);CHKERRQ(ierr); 296 297 /* store column indices (zero start index) */ 298 jj = (int*)PetscMalloc((a->s_nz+1)*bs2*sizeof(int));CHKERRQ(ierr); 299 count = 0; 300 for (i=0; i<a->mbs; i++) { 301 for (j=0; j<bs; j++) { 302 for (k=a->i[i]; k<a->i[i+1]; k++) { 303 for (l=0; l<bs; l++) { 304 jj[count++] = bs*a->j[k] + l; 305 } 306 } 307 } 308 } 309 ierr = PetscBinaryWrite(fd,jj,bs2*a->s_nz,PETSC_INT,0);CHKERRQ(ierr); 310 ierr = PetscFree(jj);CHKERRQ(ierr); 311 312 /* store nonzero values */ 313 aa = (Scalar*)PetscMalloc((a->s_nz+1)*bs2*sizeof(Scalar));CHKERRQ(ierr); 314 count = 0; 315 for (i=0; i<a->mbs; i++) { 316 for (j=0; j<bs; j++) { 317 for (k=a->i[i]; k<a->i[i+1]; k++) { 318 for (l=0; l<bs; l++) { 319 aa[count++] = a->a[bs2*k + l*bs + j]; 320 } 321 } 322 } 323 } 324 ierr = PetscBinaryWrite(fd,aa,bs2*a->s_nz,PETSC_SCALAR,0);CHKERRQ(ierr); 325 ierr = PetscFree(aa);CHKERRQ(ierr); 326 327 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 328 if (file) { 329 fprintf(file,"-matload_block_size %d\n",a->bs); 330 } 331 PetscFunctionReturn(0); 332 } 333 334 #undef __FUNC__ 335 #define __FUNC__ "MatView_SeqSBAIJ_ASCII" 336 static int MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 337 { 338 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 339 int ierr,i,j,format,bs = a->bs,k,l,bs2=a->bs2; 340 char *outputname; 341 342 PetscFunctionBegin; 343 ierr = PetscViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 344 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 345 if (format == PETSC_VIEWER_FORMAT_ASCII_INFO || format == PETSC_VIEWER_FORMAT_ASCII_INFO_LONG) { 346 ierr = PetscViewerASCIIPrintf(viewer," block size is %d\n",bs);CHKERRQ(ierr); 347 } else if (format == PETSC_VIEWER_FORMAT_ASCII_MATLAB) { 348 SETERRQ(PETSC_ERR_SUP,"Matlab format not supported"); 349 } else if (format == PETSC_VIEWER_FORMAT_ASCII_COMMON) { 350 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 351 for (i=0; i<a->mbs; i++) { 352 for (j=0; j<bs; j++) { 353 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 354 for (k=a->i[i]; k<a->i[i+1]; k++) { 355 for (l=0; l<bs; l++) { 356 #if defined(PETSC_USE_COMPLEX) 357 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 358 ierr = PetscViewerASCIIPrintf(viewer," %d %g + %g i",bs*a->j[k]+l, 359 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 360 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 361 ierr = PetscViewerASCIIPrintf(viewer," %d %g - %g i",bs*a->j[k]+l, 362 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 363 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 364 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 365 } 366 #else 367 if (a->a[bs2*k + l*bs + j] != 0.0) { 368 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 369 } 370 #endif 371 } 372 } 373 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 374 } 375 } 376 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 377 } else { 378 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 379 for (i=0; i<a->mbs; i++) { 380 for (j=0; j<bs; j++) { 381 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 382 for (k=a->i[i]; k<a->i[i+1]; k++) { 383 for (l=0; l<bs; l++) { 384 #if defined(PETSC_USE_COMPLEX) 385 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 386 ierr = PetscViewerASCIIPrintf(viewer," %d %g + %g i",bs*a->j[k]+l, 387 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 388 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 389 ierr = PetscViewerASCIIPrintf(viewer," %d %g - %g i",bs*a->j[k]+l, 390 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 391 } else { 392 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 393 } 394 #else 395 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 396 #endif 397 } 398 } 399 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 400 } 401 } 402 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 403 } 404 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 405 PetscFunctionReturn(0); 406 } 407 408 #undef __FUNC__ 409 #define __FUNC__ "MatView_SeqSBAIJ_Draw_Zoom" 410 static int MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 411 { 412 Mat A = (Mat) Aa; 413 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 414 int row,ierr,i,j,k,l,mbs=a->mbs,color,bs=a->bs,bs2=a->bs2,rank; 415 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 416 MatScalar *aa; 417 MPI_Comm comm; 418 PetscViewer viewer; 419 420 PetscFunctionBegin; 421 /* 422 This is nasty. If this is called from an originally parallel matrix 423 then all processes call this,but only the first has the matrix so the 424 rest should return immediately. 425 */ 426 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 427 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 428 if (rank) PetscFunctionReturn(0); 429 430 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 431 432 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 433 PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric"); 434 435 /* loop over matrix elements drawing boxes */ 436 color = PETSC_DRAW_BLUE; 437 for (i=0,row=0; i<mbs; i++,row+=bs) { 438 for (j=a->i[i]; j<a->i[i+1]; j++) { 439 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 440 x_l = a->j[j]*bs; x_r = x_l + 1.0; 441 aa = a->a + j*bs2; 442 for (k=0; k<bs; k++) { 443 for (l=0; l<bs; l++) { 444 if (PetscRealPart(*aa++) >= 0.) continue; 445 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 446 } 447 } 448 } 449 } 450 color = PETSC_DRAW_CYAN; 451 for (i=0,row=0; i<mbs; i++,row+=bs) { 452 for (j=a->i[i]; j<a->i[i+1]; j++) { 453 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 454 x_l = a->j[j]*bs; x_r = x_l + 1.0; 455 aa = a->a + j*bs2; 456 for (k=0; k<bs; k++) { 457 for (l=0; l<bs; l++) { 458 if (PetscRealPart(*aa++) != 0.) continue; 459 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 460 } 461 } 462 } 463 } 464 465 color = PETSC_DRAW_RED; 466 for (i=0,row=0; i<mbs; i++,row+=bs) { 467 for (j=a->i[i]; j<a->i[i+1]; j++) { 468 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 469 x_l = a->j[j]*bs; x_r = x_l + 1.0; 470 aa = a->a + j*bs2; 471 for (k=0; k<bs; k++) { 472 for (l=0; l<bs; l++) { 473 if (PetscRealPart(*aa++) <= 0.) continue; 474 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 475 } 476 } 477 } 478 } 479 PetscFunctionReturn(0); 480 } 481 482 #undef __FUNC__ 483 #define __FUNC__ "MatView_SeqSBAIJ_Draw" 484 static int MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 485 { 486 int ierr; 487 PetscReal xl,yl,xr,yr,w,h; 488 PetscDraw draw; 489 PetscTruth isnull; 490 491 PetscFunctionBegin; 492 493 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 494 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 495 496 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 497 xr = A->m; yr = A->m; h = yr/10.0; w = xr/10.0; 498 xr += w; yr += h; xl = -w; yl = -h; 499 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 500 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 501 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 502 PetscFunctionReturn(0); 503 } 504 505 #undef __FUNC__ 506 #define __FUNC__ "MatView_SeqSBAIJ" 507 int MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 508 { 509 int ierr; 510 PetscTruth issocket,isascii,isbinary,isdraw; 511 512 PetscFunctionBegin; 513 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 514 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 515 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_BINARY_VIEWER,&isbinary);CHKERRQ(ierr); 516 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_DRAW_VIEWER,&isdraw);CHKERRQ(ierr); 517 if (issocket) { 518 SETERRQ(PETSC_ERR_SUP,"Socket viewer not supported"); 519 } else if (isascii){ 520 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 521 } else if (isbinary) { 522 ierr = MatView_SeqSBAIJ_Binary(A,viewer);CHKERRQ(ierr); 523 } else if (isdraw) { 524 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 525 } else { 526 SETERRQ1(1,"Viewer type %s not supported by SeqBAIJ matrices",((PetscObject)viewer)->type_name); 527 } 528 PetscFunctionReturn(0); 529 } 530 531 532 #undef __FUNC__ 533 #define __FUNC__ "MatGetValues_SeqSBAIJ" 534 int MatGetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 535 { 536 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 537 int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 538 int *ai = a->i,*ailen = a->ilen; 539 int brow,bcol,ridx,cidx,bs=a->bs,bs2=a->bs2; 540 MatScalar *ap,*aa = a->a,zero = 0.0; 541 542 PetscFunctionBegin; 543 for (k=0; k<m; k++) { /* loop over rows */ 544 row = im[k]; brow = row/bs; 545 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative row"); 546 if (row >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 547 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 548 nrow = ailen[brow]; 549 for (l=0; l<n; l++) { /* loop over columns */ 550 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative column"); 551 if (in[l] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 552 col = in[l] ; 553 bcol = col/bs; 554 cidx = col%bs; 555 ridx = row%bs; 556 high = nrow; 557 low = 0; /* assume unsorted */ 558 while (high-low > 5) { 559 t = (low+high)/2; 560 if (rp[t] > bcol) high = t; 561 else low = t; 562 } 563 for (i=low; i<high; i++) { 564 if (rp[i] > bcol) break; 565 if (rp[i] == bcol) { 566 *v++ = ap[bs2*i+bs*cidx+ridx]; 567 goto finished; 568 } 569 } 570 *v++ = zero; 571 finished:; 572 } 573 } 574 PetscFunctionReturn(0); 575 } 576 577 578 #undef __FUNC__ 579 #define __FUNC__ "MatSetValuesBlocked_SeqSBAIJ" 580 int MatSetValuesBlocked_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 581 { 582 PetscFunctionBegin; 583 SETERRQ(1,"Function not yet written for SBAIJ format"); 584 /* PetscFunctionReturn(0); */ 585 } 586 587 #undef __FUNC__ 588 #define __FUNC__ "MatAssemblyEnd_SeqSBAIJ" 589 int MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 590 { 591 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 592 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 593 int m = A->m,*ip,N,*ailen = a->ilen; 594 int mbs = a->mbs,bs2 = a->bs2,rmax = 0,ierr; 595 MatScalar *aa = a->a,*ap; 596 597 PetscFunctionBegin; 598 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 599 600 if (m) rmax = ailen[0]; 601 for (i=1; i<mbs; i++) { 602 /* move each row back by the amount of empty slots (fshift) before it*/ 603 fshift += imax[i-1] - ailen[i-1]; 604 rmax = PetscMax(rmax,ailen[i]); 605 if (fshift) { 606 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 607 N = ailen[i]; 608 for (j=0; j<N; j++) { 609 ip[j-fshift] = ip[j]; 610 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 611 } 612 } 613 ai[i] = ai[i-1] + ailen[i-1]; 614 } 615 if (mbs) { 616 fshift += imax[mbs-1] - ailen[mbs-1]; 617 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 618 } 619 /* reset ilen and imax for each row */ 620 for (i=0; i<mbs; i++) { 621 ailen[i] = imax[i] = ai[i+1] - ai[i]; 622 } 623 a->s_nz = ai[mbs]; 624 625 /* diagonals may have moved, so kill the diagonal pointers */ 626 if (fshift && a->diag) { 627 ierr = PetscFree(a->diag);CHKERRQ(ierr); 628 PetscLogObjectMemory(A,-(m+1)*sizeof(int)); 629 a->diag = 0; 630 } 631 PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Matrix size: %d X %d, block size %d; storage space: %d unneeded, %d used\n", 632 m,A->m,a->bs,fshift*bs2,a->s_nz*bs2); 633 PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Number of mallocs during MatSetValues is %d\n", 634 a->reallocs); 635 PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Most nonzeros blocks in any row is %d\n",rmax); 636 a->reallocs = 0; 637 A->info.nz_unneeded = (PetscReal)fshift*bs2; 638 639 PetscFunctionReturn(0); 640 } 641 642 /* 643 This function returns an array of flags which indicate the locations of contiguous 644 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 645 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 646 Assume: sizes should be long enough to hold all the values. 647 */ 648 #undef __FUNC__ 649 #define __FUNC__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 650 static int MatZeroRows_SeqSBAIJ_Check_Blocks(int idx[],int n,int bs,int sizes[], int *bs_max) 651 { 652 int i,j,k,row; 653 PetscTruth flg; 654 655 PetscFunctionBegin; 656 for (i=0,j=0; i<n; j++) { 657 row = idx[i]; 658 if (row%bs!=0) { /* Not the begining of a block */ 659 sizes[j] = 1; 660 i++; 661 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 662 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 663 i++; 664 } else { /* Begining of the block, so check if the complete block exists */ 665 flg = PETSC_TRUE; 666 for (k=1; k<bs; k++) { 667 if (row+k != idx[i+k]) { /* break in the block */ 668 flg = PETSC_FALSE; 669 break; 670 } 671 } 672 if (flg == PETSC_TRUE) { /* No break in the bs */ 673 sizes[j] = bs; 674 i+= bs; 675 } else { 676 sizes[j] = 1; 677 i++; 678 } 679 } 680 } 681 *bs_max = j; 682 PetscFunctionReturn(0); 683 } 684 685 #undef __FUNC__ 686 #define __FUNC__ "MatZeroRows_SeqSBAIJ" 687 int MatZeroRows_SeqSBAIJ(Mat A,IS is,Scalar *diag) 688 { 689 Mat_SeqSBAIJ *sbaij=(Mat_SeqSBAIJ*)A->data; 690 int ierr,i,j,k,count,is_n,*is_idx,*rows; 691 int bs=sbaij->bs,bs2=sbaij->bs2,*sizes,row,bs_max; 692 Scalar zero = 0.0; 693 MatScalar *aa; 694 695 PetscFunctionBegin; 696 /* Make a copy of the IS and sort it */ 697 ierr = ISGetSize(is,&is_n);CHKERRQ(ierr); 698 ierr = ISGetIndices(is,&is_idx);CHKERRQ(ierr); 699 700 /* allocate memory for rows,sizes */ 701 ierr = PetscMalloc((3*is_n+1)*sizeof(int),& rows );CHKERRQ(ierr); 702 sizes = rows + is_n; 703 704 /* initialize copy IS values to rows, and sort them */ 705 for (i=0; i<is_n; i++) { rows[i] = is_idx[i]; } 706 ierr = PetscSortInt(is_n,rows);CHKERRQ(ierr); 707 if (sbaij->keepzeroedrows) { /* do not change nonzero structure */ 708 for (i=0; i<is_n; i++) { sizes[i] = 1; } /* sizes: size of blocks, = 1 or bs */ 709 bs_max = is_n; /* bs_max: num. of contiguous block row in the row */ 710 } else { 711 ierr = MatZeroRows_SeqSBAIJ_Check_Blocks(rows,is_n,bs,sizes,&bs_max);CHKERRQ(ierr); 712 } 713 ierr = ISRestoreIndices(is,&is_idx);CHKERRQ(ierr); 714 715 for (i=0,j=0; i<bs_max; j+=sizes[i],i++) { 716 row = rows[j]; /* row to be zeroed */ 717 if (row < 0 || row > A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %d out of range",row); 718 count = (sbaij->i[row/bs +1] - sbaij->i[row/bs])*bs; /* num. of elements in the row */ 719 aa = sbaij->a + sbaij->i[row/bs]*bs2 + (row%bs); 720 if (sizes[i] == bs && !sbaij->keepzeroedrows) { 721 if (diag) { 722 if (sbaij->ilen[row/bs] > 0) { 723 sbaij->ilen[row/bs] = 1; 724 sbaij->j[sbaij->i[row/bs]] = row/bs; 725 ierr = PetscMemzero(aa,count*bs*sizeof(MatScalar));CHKERRQ(ierr); 726 } 727 /* Now insert all the diagoanl values for this bs */ 728 for (k=0; k<bs; k++) { 729 ierr = (*A->ops->setvalues)(A,1,rows+j+k,1,rows+j+k,diag,INSERT_VALUES);CHKERRQ(ierr); 730 } 731 } else { /* (!diag) */ 732 sbaij->ilen[row/bs] = 0; 733 } /* end (!diag) */ 734 } else { /* (sizes[i] != bs), broken block */ 735 #if defined (PETSC_USE_DEBUG) 736 if (sizes[i] != 1) SETERRQ(1,"Internal Error. Value should be 1"); 737 #endif 738 for (k=0; k<count; k++) { 739 aa[0] = zero; 740 aa+=bs; 741 } 742 if (diag) { 743 ierr = (*A->ops->setvalues)(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); 744 } 745 } 746 } 747 748 ierr = PetscFree(rows);CHKERRQ(ierr); 749 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 750 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 751 PetscFunctionReturn(0); 752 } 753 754 /* Only add/insert a(i,j) with i<=j (blocks). 755 Any a(i,j) with i>j input by user is ingored. 756 */ 757 758 #undef __FUNC__ 759 #define __FUNC__ "MatSetValues_SeqSBAIJ" 760 int MatSetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 761 { 762 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 763 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 764 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 765 int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol; 766 int ridx,cidx,bs2=a->bs2,ierr; 767 MatScalar *ap,value,*aa=a->a,*bap; 768 769 PetscFunctionBegin; 770 771 for (k=0; k<m; k++) { /* loop over added rows */ 772 row = im[k]; /* row number */ 773 brow = row/bs; /* block row number */ 774 if (row < 0) continue; 775 #if defined(PETSC_USE_BOPT_g) 776 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",row,A->m); 777 #endif 778 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 779 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 780 rmax = imax[brow]; /* maximum space allocated for this row */ 781 nrow = ailen[brow]; /* actual length of this row */ 782 low = 0; 783 784 for (l=0; l<n; l++) { /* loop over added columns */ 785 if (in[l] < 0) continue; 786 #if defined(PETSC_USE_BOPT_g) 787 if (in[l] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",in[l],A->m); 788 #endif 789 col = in[l]; 790 bcol = col/bs; /* block col number */ 791 792 if (brow <= bcol){ 793 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 794 if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ 795 /* element value a(k,l) */ 796 if (roworiented) { 797 value = v[l + k*n]; 798 } else { 799 value = v[k + l*m]; 800 } 801 802 /* move pointer bap to a(k,l) quickly and add/insert value */ 803 if (!sorted) low = 0; high = nrow; 804 while (high-low > 7) { 805 t = (low+high)/2; 806 if (rp[t] > bcol) high = t; 807 else low = t; 808 } 809 for (i=low; i<high; i++) { 810 /* printf("The loop of i=low.., rp[%d]=%d\n",i,rp[i]); */ 811 if (rp[i] > bcol) break; 812 if (rp[i] == bcol) { 813 bap = ap + bs2*i + bs*cidx + ridx; 814 if (is == ADD_VALUES) *bap += value; 815 else *bap = value; 816 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 817 if (brow == bcol && ridx < cidx){ 818 bap = ap + bs2*i + bs*ridx + cidx; 819 if (is == ADD_VALUES) *bap += value; 820 else *bap = value; 821 } 822 goto noinsert1; 823 } 824 } 825 826 if (nonew == 1) goto noinsert1; 827 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 828 if (nrow >= rmax) { 829 /* there is no extra room in row, therefore enlarge */ 830 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 831 MatScalar *new_a; 832 833 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 834 835 /* Malloc new storage space */ 836 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 837 ierr = PetscMalloc(len,&( new_a ));CHKERRQ(ierr); 838 new_j = (int*)(new_a + bs2*new_nz); 839 new_i = new_j + new_nz; 840 841 /* copy over old data into new slots */ 842 for (ii=0; ii<brow+1; ii++) {new_i[ii] = ai[ii];} 843 for (ii=brow+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 844 ierr = PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int));CHKERRQ(ierr); 845 len = (new_nz - CHUNKSIZE - ai[brow] - nrow); 846 ierr = PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow,len*sizeof(int));CHKERRQ(ierr); 847 ierr = PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(MatScalar));CHKERRQ(ierr); 848 ierr = PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));CHKERRQ(ierr); 849 ierr = PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE),aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(MatScalar));CHKERRQ(ierr); 850 /* free up old matrix storage */ 851 ierr = PetscFree(a->a);CHKERRQ(ierr); 852 if (!a->singlemalloc) { 853 ierr = PetscFree(a->i);CHKERRQ(ierr); 854 ierr = PetscFree(a->j);CHKERRQ(ierr); 855 } 856 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 857 a->singlemalloc = PETSC_TRUE; 858 859 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 860 rmax = imax[brow] = imax[brow] + CHUNKSIZE; 861 PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 862 a->s_maxnz += bs2*CHUNKSIZE; 863 a->reallocs++; 864 a->s_nz++; 865 } 866 867 N = nrow++ - 1; 868 /* shift up all the later entries in this row */ 869 for (ii=N; ii>=i; ii--) { 870 rp[ii+1] = rp[ii]; 871 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 872 } 873 if (N>=i) { 874 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 875 } 876 rp[i] = bcol; 877 ap[bs2*i + bs*cidx + ridx] = value; 878 noinsert1:; 879 low = i; 880 /* } */ 881 } 882 } /* end of if .. if.. */ 883 } /* end of loop over added columns */ 884 ailen[brow] = nrow; 885 } /* end of loop over added rows */ 886 887 PetscFunctionReturn(0); 888 } 889 890 extern int MatCholeskyFactorSymbolic_SeqSBAIJ(Mat,IS,PetscReal,Mat*); 891 extern int MatCholeskyFactor_SeqSBAIJ(Mat,IS,PetscReal); 892 extern int MatIncreaseOverlap_SeqSBAIJ(Mat,int,IS*,int); 893 extern int MatGetSubMatrix_SeqSBAIJ(Mat,IS,IS,int,MatReuse,Mat*); 894 extern int MatGetSubMatrices_SeqSBAIJ(Mat,int,IS*,IS*,MatReuse,Mat**); 895 extern int MatMultTranspose_SeqSBAIJ(Mat,Vec,Vec); 896 extern int MatMultTransposeAdd_SeqSBAIJ(Mat,Vec,Vec,Vec); 897 extern int MatScale_SeqSBAIJ(Scalar*,Mat); 898 extern int MatNorm_SeqSBAIJ(Mat,NormType,PetscReal *); 899 extern int MatEqual_SeqSBAIJ(Mat,Mat,PetscTruth*); 900 extern int MatGetDiagonal_SeqSBAIJ(Mat,Vec); 901 extern int MatDiagonalScale_SeqSBAIJ(Mat,Vec,Vec); 902 extern int MatGetInfo_SeqSBAIJ(Mat,MatInfoType,MatInfo *); 903 extern int MatZeroEntries_SeqSBAIJ(Mat); 904 extern int MatGetRowMax_SeqSBAIJ(Mat,Vec); 905 906 extern int MatSolve_SeqSBAIJ_N(Mat,Vec,Vec); 907 extern int MatSolve_SeqSBAIJ_1(Mat,Vec,Vec); 908 extern int MatSolve_SeqSBAIJ_2(Mat,Vec,Vec); 909 extern int MatSolve_SeqSBAIJ_3(Mat,Vec,Vec); 910 extern int MatSolve_SeqSBAIJ_4(Mat,Vec,Vec); 911 extern int MatSolve_SeqSBAIJ_5(Mat,Vec,Vec); 912 extern int MatSolve_SeqSBAIJ_6(Mat,Vec,Vec); 913 extern int MatSolve_SeqSBAIJ_7(Mat,Vec,Vec); 914 extern int MatSolveTranspose_SeqSBAIJ_7(Mat,Vec,Vec); 915 extern int MatSolveTranspose_SeqSBAIJ_6(Mat,Vec,Vec); 916 extern int MatSolveTranspose_SeqSBAIJ_5(Mat,Vec,Vec); 917 extern int MatSolveTranspose_SeqSBAIJ_4(Mat,Vec,Vec); 918 extern int MatSolveTranspose_SeqSBAIJ_3(Mat,Vec,Vec); 919 extern int MatSolveTranspose_SeqSBAIJ_2(Mat,Vec,Vec); 920 extern int MatSolveTranspose_SeqSBAIJ_1(Mat,Vec,Vec); 921 922 extern int MatSolve_SeqSBAIJ_1_NaturalOrdering(Mat,Vec,Vec); 923 extern int MatSolve_SeqSBAIJ_2_NaturalOrdering(Mat,Vec,Vec); 924 extern int MatSolve_SeqSBAIJ_3_NaturalOrdering(Mat,Vec,Vec); 925 extern int MatSolve_SeqSBAIJ_4_NaturalOrdering(Mat,Vec,Vec); 926 extern int MatSolve_SeqSBAIJ_5_NaturalOrdering(Mat,Vec,Vec); 927 extern int MatSolve_SeqSBAIJ_6_NaturalOrdering(Mat,Vec,Vec); 928 extern int MatSolve_SeqSBAIJ_7_NaturalOrdering(Mat,Vec,Vec); 929 extern int MatSolve_SeqSBAIJ_N_NaturalOrdering(Mat,Vec,Vec); 930 931 extern int MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat,Mat*); 932 extern int MatCholeskyFactorNumeric_SeqSBAIJ_1(Mat,Mat*); 933 extern int MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat,Mat*); 934 extern int MatCholeskyFactorNumeric_SeqSBAIJ_3(Mat,Mat*); 935 extern int MatCholeskyFactorNumeric_SeqSBAIJ_4(Mat,Mat*); 936 extern int MatCholeskyFactorNumeric_SeqSBAIJ_5(Mat,Mat*); 937 extern int MatCholeskyFactorNumeric_SeqSBAIJ_6(Mat,Mat*); 938 extern int MatCholeskyFactorNumeric_SeqSBAIJ_7(Mat,Mat*); 939 940 extern int MatMult_SeqSBAIJ_1(Mat,Vec,Vec); 941 extern int MatMult_SeqSBAIJ_2(Mat,Vec,Vec); 942 extern int MatMult_SeqSBAIJ_3(Mat,Vec,Vec); 943 extern int MatMult_SeqSBAIJ_4(Mat,Vec,Vec); 944 extern int MatMult_SeqSBAIJ_5(Mat,Vec,Vec); 945 extern int MatMult_SeqSBAIJ_6(Mat,Vec,Vec); 946 extern int MatMult_SeqSBAIJ_7(Mat,Vec,Vec); 947 extern int MatMult_SeqSBAIJ_N(Mat,Vec,Vec); 948 949 extern int MatMultAdd_SeqSBAIJ_1(Mat,Vec,Vec,Vec); 950 extern int MatMultAdd_SeqSBAIJ_2(Mat,Vec,Vec,Vec); 951 extern int MatMultAdd_SeqSBAIJ_3(Mat,Vec,Vec,Vec); 952 extern int MatMultAdd_SeqSBAIJ_4(Mat,Vec,Vec,Vec); 953 extern int MatMultAdd_SeqSBAIJ_5(Mat,Vec,Vec,Vec); 954 extern int MatMultAdd_SeqSBAIJ_6(Mat,Vec,Vec,Vec); 955 extern int MatMultAdd_SeqSBAIJ_7(Mat,Vec,Vec,Vec); 956 extern int MatMultAdd_SeqSBAIJ_N(Mat,Vec,Vec,Vec); 957 958 #ifdef HAVE_MatIncompleteCholeskyFactor 959 /* modefied from MatILUFactor_SeqSBAIJ, needs further work! */ 960 #undef __FUNC__ 961 #define __FUNC__ "MatIncompleteCholeskyFactor_SeqSBAIJ" 962 int MatIncompleteCholeskyFactor_SeqSBAIJ(Mat inA,IS row,PetscReal fill,int level) 963 { 964 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 965 Mat outA; 966 int ierr; 967 PetscTruth row_identity,col_identity; 968 969 PetscFunctionBegin; 970 /* 971 if (level != 0) SETERRQ(PETSC_ERR_SUP,"Only levels = 0 supported for in-place ILU"); 972 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 973 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 974 if (!row_identity || !col_identity) { 975 SETERRQ(1,"Row and column permutations must be identity for in-place ICC"); 976 } 977 */ 978 979 outA = inA; 980 inA->factor = FACTOR_CHOLESKY; 981 982 if (!a->diag) { 983 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 984 } 985 /* 986 Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver 987 for ILU(0) factorization with natural ordering 988 */ 989 switch (a->bs) { 990 case 1: 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 a->solve_work = (Scalar*)PetscMalloc((A->m+a->bs)*sizeof(Scalar));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*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 < -2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than -2: value %d",nz); 1375 if (nnz) { 1376 for (i=0; i<mbs; i++) { 1377 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 1378 } 1379 } 1380 1381 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1382 if (!flg) { 1383 switch (bs) { 1384 case 1: 1385 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1386 B->ops->solve = MatSolve_SeqSBAIJ_1; 1387 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_1; 1388 B->ops->mult = MatMult_SeqSBAIJ_1; 1389 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1390 break; 1391 case 2: 1392 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1393 B->ops->solve = MatSolve_SeqSBAIJ_2; 1394 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_2; 1395 B->ops->mult = MatMult_SeqSBAIJ_2; 1396 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1397 break; 1398 case 3: 1399 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1400 B->ops->solve = MatSolve_SeqSBAIJ_3; 1401 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_3; 1402 B->ops->mult = MatMult_SeqSBAIJ_3; 1403 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1404 break; 1405 case 4: 1406 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1407 B->ops->solve = MatSolve_SeqSBAIJ_4; 1408 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_4; 1409 B->ops->mult = MatMult_SeqSBAIJ_4; 1410 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1411 break; 1412 case 5: 1413 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1414 B->ops->solve = MatSolve_SeqSBAIJ_5; 1415 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_5; 1416 B->ops->mult = MatMult_SeqSBAIJ_5; 1417 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1418 break; 1419 case 6: 1420 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1421 B->ops->solve = MatSolve_SeqSBAIJ_6; 1422 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_6; 1423 B->ops->mult = MatMult_SeqSBAIJ_6; 1424 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1425 break; 1426 case 7: 1427 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1428 B->ops->solve = MatSolve_SeqSBAIJ_7; 1429 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_7; 1430 B->ops->mult = MatMult_SeqSBAIJ_7; 1431 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1432 break; 1433 } 1434 } 1435 1436 b->mbs = mbs; 1437 b->nbs = mbs; 1438 ierr = PetscMalloc((mbs+1)*sizeof(int),& b->imax );CHKERRQ(ierr); 1439 if (!nnz) { 1440 if (nz == PETSC_DEFAULT) nz = 5; 1441 else if (nz <= 0) nz = 1; 1442 for (i=0; i<mbs; i++) { 1443 b->imax[i] = nz; 1444 } 1445 nz = nz*mbs; /* total nz */ 1446 } else { 1447 nz = 0; 1448 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1449 } 1450 /* s_nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1451 s_nz = nz; 1452 1453 /* allocate the matrix space */ 1454 len = s_nz*sizeof(int) + s_nz*bs2*sizeof(MatScalar) + (B->m+1)*sizeof(int); 1455 ierr = PetscMalloc(len,&( b->a ));CHKERRQ(ierr); 1456 ierr = PetscMemzero(b->a,s_nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1457 b->j = (int*)(b->a + s_nz*bs2); 1458 ierr = PetscMemzero(b->j,s_nz*sizeof(int));CHKERRQ(ierr); 1459 b->i = b->j + s_nz; 1460 b->singlemalloc = PETSC_TRUE; 1461 1462 /* pointer to beginning of each row */ 1463 b->i[0] = 0; 1464 for (i=1; i<mbs+1; i++) { 1465 b->i[i] = b->i[i-1] + b->imax[i-1]; 1466 } 1467 1468 /* b->ilen will count nonzeros in each block row so far. */ 1469 ierr = PetscMalloc((mbs+1)*sizeof(int),& b->ilen ); 1470 PetscLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1471 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1472 1473 b->bs = bs; 1474 b->bs2 = bs2; 1475 b->s_nz = 0; 1476 b->s_maxnz = s_nz*bs2; 1477 1478 b->inew = 0; 1479 b->jnew = 0; 1480 b->anew = 0; 1481 b->a2anew = 0; 1482 b->permute = PETSC_FALSE; 1483 PetscFunctionReturn(0); 1484 } 1485 1486 1487 #undef __FUNC__ 1488 #define __FUNC__ "MatCreateSeqSBAIJ" 1489 /*@C 1490 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1491 compressed row) format. For good matrix assembly performance the 1492 user should preallocate the matrix storage by setting the parameter nz 1493 (or the array nnz). By setting these parameters accurately, performance 1494 during matrix assembly can be increased by more than a factor of 50. 1495 1496 Collective on MPI_Comm 1497 1498 Input Parameters: 1499 + comm - MPI communicator, set to PETSC_COMM_SELF 1500 . bs - size of block 1501 . m - number of rows, or number of columns 1502 . nz - number of block nonzeros per block row (same for all rows) 1503 - nnz - array containing the number of block nonzeros in the various block rows 1504 (possibly different for each block row) or PETSC_NULL 1505 1506 Output Parameter: 1507 . A - the symmetric matrix 1508 1509 Options Database Keys: 1510 . -mat_no_unroll - uses code that does not unroll the loops in the 1511 block calculations (much slower) 1512 . -mat_block_size - size of the blocks to use 1513 1514 Level: intermediate 1515 1516 Notes: 1517 The block AIJ format is compatible with standard Fortran 77 1518 storage. That is, the stored row and column indices can begin at 1519 either one (as in Fortran) or zero. See the users' manual for details. 1520 1521 Specify the preallocated storage with either nz or nnz (not both). 1522 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1523 allocation. For additional details, see the users manual chapter on 1524 matrices. 1525 1526 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1527 @*/ 1528 int MatCreateSeqSBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1529 { 1530 int ierr; 1531 1532 PetscFunctionBegin; 1533 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1534 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1535 ierr = MatSeqSBAIJSetPreallocation(*A,bs,nz,nnz);CHKERRQ(ierr); 1536 PetscFunctionReturn(0); 1537 } 1538 1539 #undef __FUNC__ 1540 #define __FUNC__ "MatDuplicate_SeqSBAIJ" 1541 int MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1542 { 1543 Mat C; 1544 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1545 int i,len,mbs = a->mbs,nz = a->s_nz,bs2 =a->bs2,ierr; 1546 1547 PetscFunctionBegin; 1548 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1549 1550 *B = 0; 1551 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 1552 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 1553 c = (Mat_SeqSBAIJ*)C->data; 1554 1555 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1556 C->preallocated = PETSC_TRUE; 1557 C->factor = A->factor; 1558 c->row = 0; 1559 c->icol = 0; 1560 c->saved_values = 0; 1561 c->keepzeroedrows = a->keepzeroedrows; 1562 C->assembled = PETSC_TRUE; 1563 1564 c->bs = a->bs; 1565 c->bs2 = a->bs2; 1566 c->mbs = a->mbs; 1567 c->nbs = a->nbs; 1568 1569 ierr = PetscMalloc((mbs+1)*sizeof(int),& c->imax );CHKERRQ(ierr); 1570 ierr = PetscMalloc((mbs+1)*sizeof(int),& c->ilen );CHKERRQ(ierr); 1571 for (i=0; i<mbs; i++) { 1572 c->imax[i] = a->imax[i]; 1573 c->ilen[i] = a->ilen[i]; 1574 } 1575 1576 /* allocate the matrix space */ 1577 c->singlemalloc = PETSC_TRUE; 1578 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1579 ierr = PetscMalloc(len,&( c->a ));CHKERRQ(ierr); 1580 c->j = (int*)(c->a + nz*bs2); 1581 c->i = c->j + nz; 1582 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));CHKERRQ(ierr); 1583 if (mbs > 0) { 1584 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(int));CHKERRQ(ierr); 1585 if (cpvalues == MAT_COPY_VALUES) { 1586 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1587 } else { 1588 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1589 } 1590 } 1591 1592 PetscLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1593 c->sorted = a->sorted; 1594 c->roworiented = a->roworiented; 1595 c->nonew = a->nonew; 1596 1597 if (a->diag) { 1598 ierr = PetscMalloc((mbs+1)*sizeof(int),& c->diag );CHKERRQ(ierr); 1599 PetscLogObjectMemory(C,(mbs+1)*sizeof(int)); 1600 for (i=0; i<mbs; i++) { 1601 c->diag[i] = a->diag[i]; 1602 } 1603 } else c->diag = 0; 1604 c->s_nz = a->s_nz; 1605 c->s_maxnz = a->s_maxnz; 1606 c->solve_work = 0; 1607 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1608 c->mult_work = 0; 1609 *B = C; 1610 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1611 PetscFunctionReturn(0); 1612 } 1613 1614 EXTERN_C_BEGIN 1615 #undef __FUNC__ 1616 #define __FUNC__ "MatLoad_SeqSBAIJ" 1617 int MatLoad_SeqSBAIJ(PetscViewer viewer,MatType type,Mat *A) 1618 { 1619 Mat_SeqSBAIJ *a; 1620 Mat B; 1621 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1; 1622 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,*s_browlengths,maskcount; 1623 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1624 int *masked,nmask,tmp,bs2,ishift; 1625 Scalar *aa; 1626 MPI_Comm comm = ((PetscObject)viewer)->comm; 1627 1628 PetscFunctionBegin; 1629 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1630 bs2 = bs*bs; 1631 1632 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1633 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1634 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1635 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1636 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1637 M = header[1]; N = header[2]; nz = header[3]; 1638 1639 if (header[3] < 0) { 1640 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 1641 } 1642 1643 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1644 1645 /* 1646 This code adds extra rows to make sure the number of rows is 1647 divisible by the blocksize 1648 */ 1649 mbs = M/bs; 1650 extra_rows = bs - M + bs*(mbs); 1651 if (extra_rows == bs) extra_rows = 0; 1652 else mbs++; 1653 if (extra_rows) { 1654 PetscLogInfo(0,"MatLoad_SeqSBAIJ:Padding loaded matrix to match blocksize\n"); 1655 } 1656 1657 /* read in row lengths */ 1658 ierr = PetscMalloc((M+extra_rows)*sizeof(int),& rowlengths );CHKERRQ(ierr); 1659 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1660 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1661 1662 /* read in column indices */ 1663 ierr = PetscMalloc((nz+extra_rows)*sizeof(int),& jj );CHKERRQ(ierr); 1664 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1665 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1666 1667 /* loop over row lengths determining block row lengths */ 1668 ierr = PetscMalloc(mbs*sizeof(int),&( browlengths ));CHKERRQ(ierr); 1669 ierr = PetscMalloc(mbs*sizeof(int),&( s_browlengths ));CHKERRQ(ierr); 1670 ierr = PetscMemzero(s_browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1671 ierr = PetscMalloc(2*mbs*sizeof(int),&( mask ));CHKERRQ(ierr); 1672 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1673 masked = mask + mbs; 1674 rowcount = 0; nzcount = 0; 1675 for (i=0; i<mbs; i++) { 1676 nmask = 0; 1677 for (j=0; j<bs; j++) { 1678 kmax = rowlengths[rowcount]; 1679 for (k=0; k<kmax; k++) { 1680 tmp = jj[nzcount++]/bs; /* block col. index */ 1681 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 1682 } 1683 rowcount++; 1684 } 1685 s_browlengths[i] += nmask; 1686 browlengths[i] = 2*s_browlengths[i]; 1687 1688 /* zero out the mask elements we set */ 1689 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1690 } 1691 1692 /* create our matrix */ 1693 ierr = MatCreateSeqSBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1694 B = *A; 1695 a = (Mat_SeqSBAIJ*)B->data; 1696 1697 /* set matrix "i" values */ 1698 a->i[0] = 0; 1699 for (i=1; i<= mbs; i++) { 1700 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 1701 a->ilen[i-1] = s_browlengths[i-1]; 1702 } 1703 a->s_nz = 0; 1704 for (i=0; i<mbs; i++) a->s_nz += s_browlengths[i]; 1705 1706 /* read in nonzero values */ 1707 ierr = PetscMalloc((nz+extra_rows)*sizeof(Scalar),& aa );CHKERRQ(ierr); 1708 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1709 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1710 1711 /* set "a" and "j" values into matrix */ 1712 nzcount = 0; jcount = 0; 1713 for (i=0; i<mbs; i++) { 1714 nzcountb = nzcount; 1715 nmask = 0; 1716 for (j=0; j<bs; j++) { 1717 kmax = rowlengths[i*bs+j]; 1718 for (k=0; k<kmax; k++) { 1719 tmp = jj[nzcount++]/bs; /* block col. index */ 1720 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 1721 } 1722 rowcount++; 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