1 /*$Id: sbaij.c,v 1.43 2000/11/01 17:34:45 hzhang Exp hzhang $*/ 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 diag = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(diag); 47 PLogObjectMemory(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 PLogObjectState((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 PLogInfo(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 *v = (Scalar*)PetscMalloc((*ncols+row)*sizeof(Scalar));CHKPTRQ(*v); 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 *cols = (int*)PetscMalloc((*ncols+row)*sizeof(int));CHKPTRQ(*cols); 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,Viewer 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 = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 280 col_lens = (int*)PetscMalloc((4+A->m)*sizeof(int));CHKPTRQ(col_lens); 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));CHKPTRQ(jj); 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));CHKPTRQ(aa); 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 = ViewerBinaryGetInfoPointer(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,Viewer 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 = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 344 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 345 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 346 ierr = ViewerASCIIPrintf(viewer," block size is %d\n",bs);CHKERRQ(ierr); 347 } else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 348 SETERRQ(PETSC_ERR_SUP,"Matlab format not supported"); 349 } else if (format == VIEWER_FORMAT_ASCII_COMMON) { 350 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 351 for (i=0; i<a->mbs; i++) { 352 for (j=0; j<bs; j++) { 353 ierr = ViewerASCIIPrintf(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 = ViewerASCIIPrintf(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 = ViewerASCIIPrintf(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 = ViewerASCIIPrintf(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 = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 369 } 370 #endif 371 } 372 } 373 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 374 } 375 } 376 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 377 } else { 378 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 379 for (i=0; i<a->mbs; i++) { 380 for (j=0; j<bs; j++) { 381 ierr = ViewerASCIIPrintf(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 = ViewerASCIIPrintf(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 = ViewerASCIIPrintf(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 = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 393 } 394 #else 395 ierr = ViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 396 #endif 397 } 398 } 399 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 400 } 401 } 402 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 403 } 404 ierr = ViewerFlush(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(Draw 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 Viewer 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 = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 433 DrawString(draw, .3*(xl+xr), .3*(yl+yr), DRAW_BLACK, "symmetric"); 434 435 /* loop over matrix elements drawing boxes */ 436 color = 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 = DrawRectangle(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 = 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 = DrawRectangle(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 = 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 = DrawRectangle(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,Viewer viewer) 485 { 486 int ierr; 487 PetscReal xl,yl,xr,yr,w,h; 488 Draw draw; 489 PetscTruth isnull; 490 491 PetscFunctionBegin; 492 493 ierr = ViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 494 ierr = DrawIsNull(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 = DrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 500 ierr = DrawZoom(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,Viewer viewer) 508 { 509 int ierr; 510 PetscTruth issocket,isascii,isbinary,isdraw; 511 512 PetscFunctionBegin; 513 ierr = PetscTypeCompare((PetscObject)viewer,SOCKET_VIEWER,&issocket);CHKERRQ(ierr); 514 ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr); 515 ierr = PetscTypeCompare((PetscObject)viewer,BINARY_VIEWER,&isbinary);CHKERRQ(ierr); 516 ierr = PetscTypeCompare((PetscObject)viewer,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 PLogObjectMemory(A,-(m+1)*sizeof(int)); 629 a->diag = 0; 630 } 631 PLogInfo(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 PLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Number of mallocs during MatSetValues is %d\n", 634 a->reallocs); 635 PLogInfo(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 rows = (int*)PetscMalloc((3*is_n+1)*sizeof(int));CHKPTRQ(rows); 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 new_a = (MatScalar*)PetscMalloc(len);CHKPTRQ(new_a); 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 PLogObjectMemory(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 Ploginfo(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 PLogInfo(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 PLogInfo(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 PLogInfo(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 PLogInfo(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 PLogInfo(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 PLogInfo(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 PLogObjectParent(inA,a->icol); 1038 1039 if (!a->solve_work) { 1040 a->solve_work = (Scalar*)PetscMalloc((A->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1041 PLogObjectMemory(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__ /*<a name="MatSetUpPreallocation_SeqSBAIJ"></a>*/"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 aij->saved_values = (Scalar*)PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 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 B->data = (void*)(b = PetscNew(Mat_SeqSBAIJ));CHKPTRQ(b); 1276 ierr = PetscMemzero(b,sizeof(Mat_SeqSBAIJ));CHKERRQ(ierr); 1277 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1278 B->ops->destroy = MatDestroy_SeqSBAIJ; 1279 B->ops->view = MatView_SeqSBAIJ; 1280 B->factor = 0; 1281 B->lupivotthreshold = 1.0; 1282 B->mapping = 0; 1283 b->row = 0; 1284 b->icol = 0; 1285 b->reallocs = 0; 1286 b->saved_values = 0; 1287 1288 ierr = MapCreateMPI(B->comm,B->m,B->M,&B->rmap);CHKERRQ(ierr); 1289 ierr = MapCreateMPI(B->comm,B->n,B->N,&B->cmap);CHKERRQ(ierr); 1290 1291 b->sorted = PETSC_FALSE; 1292 b->roworiented = PETSC_TRUE; 1293 b->nonew = 0; 1294 b->diag = 0; 1295 b->solve_work = 0; 1296 b->mult_work = 0; 1297 b->spptr = 0; 1298 b->keepzeroedrows = PETSC_FALSE; 1299 1300 b->inew = 0; 1301 b->jnew = 0; 1302 b->anew = 0; 1303 b->a2anew = 0; 1304 b->permute = PETSC_FALSE; 1305 1306 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1307 "MatStoreValues_SeqSBAIJ", 1308 (void*)MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1309 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1310 "MatRetrieveValues_SeqSBAIJ", 1311 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1312 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1313 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1314 (void*)MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1315 PetscFunctionReturn(0); 1316 } 1317 EXTERN_C_END 1318 1319 #undef __FUNC__ 1320 #define __FUNC__ "MatSeqSBAIJSetPreallocation" 1321 /*@C 1322 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1323 compressed row) format. For good matrix assembly performance the 1324 user should preallocate the matrix storage by setting the parameter nz 1325 (or the array nnz). By setting these parameters accurately, performance 1326 during matrix assembly can be increased by more than a factor of 50. 1327 1328 Collective on Mat 1329 1330 Input Parameters: 1331 + A - the symmetric matrix 1332 . bs - size of block 1333 . nz - number of block nonzeros per block row (same for all rows) 1334 - nnz - array containing the number of block nonzeros in the various block rows 1335 (possibly different for each block row) or PETSC_NULL 1336 1337 Options Database Keys: 1338 . -mat_no_unroll - uses code that does not unroll the loops in the 1339 block calculations (much slower) 1340 . -mat_block_size - size of the blocks to use 1341 1342 Level: intermediate 1343 1344 Notes: 1345 The block AIJ format is compatible with standard Fortran 77 1346 storage. That is, the stored row and column indices can begin at 1347 either one (as in Fortran) or zero. See the users' manual for details. 1348 1349 Specify the preallocated storage with either nz or nnz (not both). 1350 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1351 allocation. For additional details, see the users manual chapter on 1352 matrices. 1353 1354 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1355 @*/ 1356 int MatSeqSBAIJSetPreallocation(Mat B,int bs,int nz,int *nnz) 1357 { 1358 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1359 int i,len,ierr,mbs,bs2; 1360 PetscTruth flg; 1361 int s_nz; 1362 1363 PetscFunctionBegin; 1364 B->preallocated = PETSC_TRUE; 1365 ierr = OptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1366 mbs = B->m/bs; 1367 bs2 = bs*bs; 1368 1369 if (mbs*bs != B->m) { 1370 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1371 } 1372 1373 if (nz < -2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than -2: value %d",nz); 1374 if (nnz) { 1375 for (i=0; i<mbs; i++) { 1376 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 1377 } 1378 } 1379 1380 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1381 if (!flg) { 1382 switch (bs) { 1383 case 1: 1384 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1385 B->ops->solve = MatSolve_SeqSBAIJ_1; 1386 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_1; 1387 B->ops->mult = MatMult_SeqSBAIJ_1; 1388 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1389 break; 1390 case 2: 1391 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1392 B->ops->solve = MatSolve_SeqSBAIJ_2; 1393 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_2; 1394 B->ops->mult = MatMult_SeqSBAIJ_2; 1395 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1396 break; 1397 case 3: 1398 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1399 B->ops->solve = MatSolve_SeqSBAIJ_3; 1400 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_3; 1401 B->ops->mult = MatMult_SeqSBAIJ_3; 1402 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1403 break; 1404 case 4: 1405 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1406 B->ops->solve = MatSolve_SeqSBAIJ_4; 1407 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_4; 1408 B->ops->mult = MatMult_SeqSBAIJ_4; 1409 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1410 break; 1411 case 5: 1412 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1413 B->ops->solve = MatSolve_SeqSBAIJ_5; 1414 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_5; 1415 B->ops->mult = MatMult_SeqSBAIJ_5; 1416 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1417 break; 1418 case 6: 1419 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1420 B->ops->solve = MatSolve_SeqSBAIJ_6; 1421 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_6; 1422 B->ops->mult = MatMult_SeqSBAIJ_6; 1423 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1424 break; 1425 case 7: 1426 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1427 B->ops->solve = MatSolve_SeqSBAIJ_7; 1428 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_7; 1429 B->ops->mult = MatMult_SeqSBAIJ_7; 1430 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1431 break; 1432 } 1433 } 1434 1435 b->mbs = mbs; 1436 b->nbs = mbs; 1437 b->imax = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(b->imax); 1438 if (!nnz) { 1439 if (nz == PETSC_DEFAULT) nz = 5; 1440 else if (nz <= 0) nz = 1; 1441 for (i=0; i<mbs; i++) { 1442 b->imax[i] = nz; 1443 } 1444 nz = nz*mbs; /* total nz */ 1445 } else { 1446 nz = 0; 1447 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1448 } 1449 /* s_nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1450 s_nz = nz; 1451 1452 /* allocate the matrix space */ 1453 len = s_nz*sizeof(int) + s_nz*bs2*sizeof(MatScalar) + (B->m+1)*sizeof(int); 1454 b->a = (MatScalar*)PetscMalloc(len);CHKPTRQ(b->a); 1455 ierr = PetscMemzero(b->a,s_nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1456 b->j = (int*)(b->a + s_nz*bs2); 1457 ierr = PetscMemzero(b->j,s_nz*sizeof(int));CHKERRQ(ierr); 1458 b->i = b->j + s_nz; 1459 b->singlemalloc = PETSC_TRUE; 1460 1461 /* pointer to beginning of each row */ 1462 b->i[0] = 0; 1463 for (i=1; i<mbs+1; i++) { 1464 b->i[i] = b->i[i-1] + b->imax[i-1]; 1465 } 1466 1467 /* b->ilen will count nonzeros in each block row so far. */ 1468 b->ilen = (int*)PetscMalloc((mbs+1)*sizeof(int)); 1469 PLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1470 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1471 1472 b->bs = bs; 1473 b->bs2 = bs2; 1474 b->s_nz = 0; 1475 b->s_maxnz = s_nz*bs2; 1476 1477 b->inew = 0; 1478 b->jnew = 0; 1479 b->anew = 0; 1480 b->a2anew = 0; 1481 b->permute = PETSC_FALSE; 1482 PetscFunctionReturn(0); 1483 } 1484 1485 1486 #undef __FUNC__ 1487 #define __FUNC__ "MatCreateSeqSBAIJ" 1488 /*@C 1489 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1490 compressed row) format. For good matrix assembly performance the 1491 user should preallocate the matrix storage by setting the parameter nz 1492 (or the array nnz). By setting these parameters accurately, performance 1493 during matrix assembly can be increased by more than a factor of 50. 1494 1495 Collective on MPI_Comm 1496 1497 Input Parameters: 1498 + comm - MPI communicator, set to PETSC_COMM_SELF 1499 . bs - size of block 1500 . m - number of rows, or number of columns 1501 . nz - number of block nonzeros per block row (same for all rows) 1502 - nnz - array containing the number of block nonzeros in the various block rows 1503 (possibly different for each block row) or PETSC_NULL 1504 1505 Output Parameter: 1506 . A - the symmetric matrix 1507 1508 Options Database Keys: 1509 . -mat_no_unroll - uses code that does not unroll the loops in the 1510 block calculations (much slower) 1511 . -mat_block_size - size of the blocks to use 1512 1513 Level: intermediate 1514 1515 Notes: 1516 The block AIJ format is compatible with standard Fortran 77 1517 storage. That is, the stored row and column indices can begin at 1518 either one (as in Fortran) or zero. See the users' manual for details. 1519 1520 Specify the preallocated storage with either nz or nnz (not both). 1521 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1522 allocation. For additional details, see the users manual chapter on 1523 matrices. 1524 1525 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1526 @*/ 1527 int MatCreateSeqSBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1528 { 1529 int ierr; 1530 1531 PetscFunctionBegin; 1532 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1533 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1534 ierr = MatSeqSBAIJSetPreallocation(*A,bs,nz,nnz);CHKERRQ(ierr); 1535 PetscFunctionReturn(0); 1536 } 1537 1538 #undef __FUNC__ 1539 #define __FUNC__ "MatDuplicate_SeqSBAIJ" 1540 int MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1541 { 1542 Mat C; 1543 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1544 int i,len,mbs = a->mbs,nz = a->s_nz,bs2 =a->bs2,ierr; 1545 1546 PetscFunctionBegin; 1547 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1548 1549 *B = 0; 1550 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 1551 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 1552 c = (Mat_SeqSBAIJ*)C->data; 1553 1554 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1555 C->preallocated = PETSC_TRUE; 1556 C->factor = A->factor; 1557 c->row = 0; 1558 c->icol = 0; 1559 c->saved_values = 0; 1560 c->keepzeroedrows = a->keepzeroedrows; 1561 C->assembled = PETSC_TRUE; 1562 1563 c->bs = a->bs; 1564 c->bs2 = a->bs2; 1565 c->mbs = a->mbs; 1566 c->nbs = a->nbs; 1567 1568 c->imax = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->imax); 1569 c->ilen = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->ilen); 1570 for (i=0; i<mbs; i++) { 1571 c->imax[i] = a->imax[i]; 1572 c->ilen[i] = a->ilen[i]; 1573 } 1574 1575 /* allocate the matrix space */ 1576 c->singlemalloc = PETSC_TRUE; 1577 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1578 c->a = (MatScalar*)PetscMalloc(len);CHKPTRQ(c->a); 1579 c->j = (int*)(c->a + nz*bs2); 1580 c->i = c->j + nz; 1581 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));CHKERRQ(ierr); 1582 if (mbs > 0) { 1583 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(int));CHKERRQ(ierr); 1584 if (cpvalues == MAT_COPY_VALUES) { 1585 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1586 } else { 1587 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1588 } 1589 } 1590 1591 PLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1592 c->sorted = a->sorted; 1593 c->roworiented = a->roworiented; 1594 c->nonew = a->nonew; 1595 1596 if (a->diag) { 1597 c->diag = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->diag); 1598 PLogObjectMemory(C,(mbs+1)*sizeof(int)); 1599 for (i=0; i<mbs; i++) { 1600 c->diag[i] = a->diag[i]; 1601 } 1602 } else c->diag = 0; 1603 c->s_nz = a->s_nz; 1604 c->s_maxnz = a->s_maxnz; 1605 c->solve_work = 0; 1606 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1607 c->mult_work = 0; 1608 *B = C; 1609 ierr = FListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1610 PetscFunctionReturn(0); 1611 } 1612 1613 EXTERN_C_BEGIN 1614 #undef __FUNC__ 1615 #define __FUNC__ "MatLoad_SeqSBAIJ" 1616 int MatLoad_SeqSBAIJ(Viewer viewer,MatType type,Mat *A) 1617 { 1618 Mat_SeqSBAIJ *a; 1619 Mat B; 1620 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1; 1621 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,*s_browlengths,maskcount; 1622 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1623 int *masked,nmask,tmp,bs2,ishift; 1624 Scalar *aa; 1625 MPI_Comm comm = ((PetscObject)viewer)->comm; 1626 1627 PetscFunctionBegin; 1628 ierr = OptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1629 bs2 = bs*bs; 1630 1631 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1632 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1633 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1634 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1635 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1636 M = header[1]; N = header[2]; nz = header[3]; 1637 1638 if (header[3] < 0) { 1639 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 1640 } 1641 1642 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1643 1644 /* 1645 This code adds extra rows to make sure the number of rows is 1646 divisible by the blocksize 1647 */ 1648 mbs = M/bs; 1649 extra_rows = bs - M + bs*(mbs); 1650 if (extra_rows == bs) extra_rows = 0; 1651 else mbs++; 1652 if (extra_rows) { 1653 PLogInfo(0,"MatLoad_SeqSBAIJ:Padding loaded matrix to match blocksize\n"); 1654 } 1655 1656 /* read in row lengths */ 1657 rowlengths = (int*)PetscMalloc((M+extra_rows)*sizeof(int));CHKPTRQ(rowlengths); 1658 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1659 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1660 1661 /* read in column indices */ 1662 jj = (int*)PetscMalloc((nz+extra_rows)*sizeof(int));CHKPTRQ(jj); 1663 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1664 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1665 1666 /* loop over row lengths determining block row lengths */ 1667 browlengths = (int*)PetscMalloc(mbs*sizeof(int));CHKPTRQ(browlengths); 1668 s_browlengths = (int*)PetscMalloc(mbs*sizeof(int));CHKPTRQ(s_browlengths); 1669 ierr = PetscMemzero(s_browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1670 mask = (int*)PetscMalloc(2*mbs*sizeof(int));CHKPTRQ(mask); 1671 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1672 masked = mask + mbs; 1673 rowcount = 0; nzcount = 0; 1674 for (i=0; i<mbs; i++) { 1675 nmask = 0; 1676 for (j=0; j<bs; j++) { 1677 kmax = rowlengths[rowcount]; 1678 for (k=0; k<kmax; k++) { 1679 tmp = jj[nzcount++]/bs; /* block col. index */ 1680 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 1681 } 1682 rowcount++; 1683 } 1684 s_browlengths[i] += nmask; 1685 browlengths[i] = 2*s_browlengths[i]; 1686 1687 /* zero out the mask elements we set */ 1688 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1689 } 1690 1691 /* create our matrix */ 1692 ierr = MatCreateSeqSBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1693 B = *A; 1694 a = (Mat_SeqSBAIJ*)B->data; 1695 1696 /* set matrix "i" values */ 1697 a->i[0] = 0; 1698 for (i=1; i<= mbs; i++) { 1699 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 1700 a->ilen[i-1] = s_browlengths[i-1]; 1701 } 1702 a->s_nz = 0; 1703 for (i=0; i<mbs; i++) a->s_nz += s_browlengths[i]; 1704 1705 /* read in nonzero values */ 1706 aa = (Scalar*)PetscMalloc((nz+extra_rows)*sizeof(Scalar));CHKPTRQ(aa); 1707 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1708 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1709 1710 /* set "a" and "j" values into matrix */ 1711 nzcount = 0; jcount = 0; 1712 for (i=0; i<mbs; i++) { 1713 nzcountb = nzcount; 1714 nmask = 0; 1715 for (j=0; j<bs; j++) { 1716 kmax = rowlengths[i*bs+j]; 1717 for (k=0; k<kmax; k++) { 1718 tmp = jj[nzcount++]/bs; /* block col. index */ 1719 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 1720 } 1721 rowcount++; 1722 } 1723 /* sort the masked values */ 1724 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1725 1726 /* set "j" values into matrix */ 1727 maskcount = 1; 1728 for (j=0; j<nmask; j++) { 1729 a->j[jcount++] = masked[j]; 1730 mask[masked[j]] = maskcount++; 1731 } 1732 1733 /* set "a" values into matrix */ 1734 ishift = bs2*a->i[i]; 1735 for (j=0; j<bs; j++) { 1736 kmax = rowlengths[i*bs+j]; 1737 for (k=0; k<kmax; k++) { 1738 tmp = jj[nzcountb]/bs ; /* block col. index */ 1739 if (tmp >= i){ 1740 block = mask[tmp] - 1; 1741 point = jj[nzcountb] - bs*tmp; 1742 idx = ishift + bs2*block + j + bs*point; 1743 a->a[idx] = aa[nzcountb]; 1744 } 1745 nzcountb++; 1746 } 1747 } 1748 /* zero out the mask elements we set */ 1749 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1750 } 1751 if (jcount != a->s_nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 1752 1753 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1754 ierr = PetscFree(browlengths);CHKERRQ(ierr); 1755 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 1756 ierr = PetscFree(aa);CHKERRQ(ierr); 1757 ierr = PetscFree(jj);CHKERRQ(ierr); 1758 ierr = PetscFree(mask);CHKERRQ(ierr); 1759 1760 B->assembled = PETSC_TRUE; 1761 ierr = MatView_Private(B);CHKERRQ(ierr); 1762 PetscFunctionReturn(0); 1763 } 1764 EXTERN_C_END 1765 1766 1767 1768