1 /*$Id: sbaij.c,v 1.41 2000/10/30 18:47:39 hzhang Exp balay $*/ 2 3 /* 4 Defines the basic matrix operations for the BAIJ (compressed row) 5 matrix storage format. 6 */ 7 #include "petscsys.h" /*I "petscmat.h" I*/ 8 #include "src/mat/impls/baij/seq/baij.h" 9 #include "src/vec/vecimpl.h" 10 #include "src/inline/spops.h" 11 #include "src/mat/impls/sbaij/seq/sbaij.h" 12 13 #define CHUNKSIZE 10 14 15 /* 16 Checks for missing diagonals 17 */ 18 #undef __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 939 extern int MatMult_SeqSBAIJ_1(Mat,Vec,Vec); 940 extern int MatMult_SeqSBAIJ_2(Mat,Vec,Vec); 941 extern int MatMult_SeqSBAIJ_3(Mat,Vec,Vec); 942 extern int MatMult_SeqSBAIJ_4(Mat,Vec,Vec); 943 extern int MatMult_SeqSBAIJ_5(Mat,Vec,Vec); 944 extern int MatMult_SeqSBAIJ_6(Mat,Vec,Vec); 945 extern int MatMult_SeqSBAIJ_7(Mat,Vec,Vec); 946 extern int MatMult_SeqSBAIJ_N(Mat,Vec,Vec); 947 948 extern int MatMultAdd_SeqSBAIJ_1(Mat,Vec,Vec,Vec); 949 extern int MatMultAdd_SeqSBAIJ_2(Mat,Vec,Vec,Vec); 950 extern int MatMultAdd_SeqSBAIJ_3(Mat,Vec,Vec,Vec); 951 extern int MatMultAdd_SeqSBAIJ_4(Mat,Vec,Vec,Vec); 952 extern int MatMultAdd_SeqSBAIJ_5(Mat,Vec,Vec,Vec); 953 extern int MatMultAdd_SeqSBAIJ_6(Mat,Vec,Vec,Vec); 954 extern int MatMultAdd_SeqSBAIJ_7(Mat,Vec,Vec,Vec); 955 extern int MatMultAdd_SeqSBAIJ_N(Mat,Vec,Vec,Vec); 956 957 #ifdef HAVE_MatIncompleteCholeskyFactor 958 /* modefied from MatILUFactor_SeqSBAIJ, needs further work! */ 959 #undef __FUNC__ 960 #define __FUNC__ "MatIncompleteCholeskyFactor_SeqSBAIJ" 961 int MatIncompleteCholeskyFactor_SeqSBAIJ(Mat inA,IS row,PetscReal fill,int level) 962 { 963 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 964 Mat outA; 965 int ierr; 966 PetscTruth row_identity,col_identity; 967 968 PetscFunctionBegin; 969 /* 970 if (level != 0) SETERRQ(PETSC_ERR_SUP,"Only levels = 0 supported for in-place ILU"); 971 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 972 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 973 if (!row_identity || !col_identity) { 974 SETERRQ(1,"Row and column permutations must be identity for in-place ICC"); 975 } 976 */ 977 978 outA = inA; 979 inA->factor = FACTOR_CHOLESKY; 980 981 if (!a->diag) { 982 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 983 } 984 /* 985 Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver 986 for ILU(0) factorization with natural ordering 987 */ 988 switch (a->bs) { 989 case 1: 990 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering; 991 Ploginfo(inA,"MatIncompleteCholeskyFactor_SeqSBAIJ:Using special in-place natural ordering solvetrans BS=1\n"); 992 case 2: 993 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 994 inA->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering; 995 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering; 996 PLogInfo(inA,"MatIncompleteCholeskyFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=2\n"); 997 break; 998 case 3: 999 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1000 inA->ops->solve = MatSolve_SeqSBAIJ_3_NaturalOrdering; 1001 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_3_NaturalOrdering; 1002 PLogInfo(inA,"MatIncompleteCholeskyFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=3\n"); 1003 break; 1004 case 4: 1005 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1006 inA->ops->solve = MatSolve_SeqSBAIJ_4_NaturalOrdering; 1007 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_4_NaturalOrdering; 1008 PLogInfo(inA,"MatIncompleteCholeskyFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=4\n"); 1009 break; 1010 case 5: 1011 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1012 inA->ops->solve = MatSolve_SeqSBAIJ_5_NaturalOrdering; 1013 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_5_NaturalOrdering; 1014 PLogInfo(inA,"MatIncompleteCholeskyFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=5\n"); 1015 break; 1016 case 6: 1017 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1018 inA->ops->solve = MatSolve_SeqSBAIJ_6_NaturalOrdering; 1019 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_6_NaturalOrdering; 1020 PLogInfo(inA,"MatIncompleteCholeskyFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=6\n"); 1021 break; 1022 case 7: 1023 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1024 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_7_NaturalOrdering; 1025 inA->ops->solve = MatSolve_SeqSBAIJ_7_NaturalOrdering; 1026 PLogInfo(inA,"MatIncompleteCholeskyFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=7\n"); 1027 break; 1028 default: 1029 a->row = row; 1030 a->icol = col; 1031 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1032 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1033 1034 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1035 ierr = ISInvertPermutation(col,PETSC_DECIDE, &(a->icol));CHKERRQ(ierr); 1036 PLogObjectParent(inA,a->icol); 1037 1038 if (!a->solve_work) { 1039 a->solve_work = (Scalar*)PetscMalloc((A->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1040 PLogObjectMemory(inA,(A->m+a->bs)*sizeof(Scalar)); 1041 } 1042 } 1043 1044 ierr = MatCholeskyFactorNumeric(inA,&outA);CHKERRQ(ierr); 1045 1046 PetscFunctionReturn(0); 1047 } 1048 #endif 1049 1050 #undef __FUNC__ 1051 #define __FUNC__ "MatPrintHelp_SeqSBAIJ" 1052 int MatPrintHelp_SeqSBAIJ(Mat A) 1053 { 1054 static PetscTruth called = PETSC_FALSE; 1055 MPI_Comm comm = A->comm; 1056 int ierr; 1057 1058 PetscFunctionBegin; 1059 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1060 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQSBAIJ and MATMPISBAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1061 ierr = (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n");CHKERRQ(ierr); 1062 PetscFunctionReturn(0); 1063 } 1064 1065 EXTERN_C_BEGIN 1066 #undef __FUNC__ 1067 #define __FUNC__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1068 int MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,int *indices) 1069 { 1070 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1071 int i,nz,n; 1072 1073 PetscFunctionBegin; 1074 nz = baij->s_maxnz; 1075 n = mat->n; 1076 for (i=0; i<nz; i++) { 1077 baij->j[i] = indices[i]; 1078 } 1079 baij->s_nz = nz; 1080 for (i=0; i<n; i++) { 1081 baij->ilen[i] = baij->imax[i]; 1082 } 1083 1084 PetscFunctionReturn(0); 1085 } 1086 EXTERN_C_END 1087 1088 #undef __FUNC__ 1089 #define __FUNC__ "MatSeqSBAIJSetColumnIndices" 1090 /*@ 1091 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1092 in the matrix. 1093 1094 Input Parameters: 1095 + mat - the SeqSBAIJ matrix 1096 - indices - the column indices 1097 1098 Level: advanced 1099 1100 Notes: 1101 This can be called if you have precomputed the nonzero structure of the 1102 matrix and want to provide it to the matrix object to improve the performance 1103 of the MatSetValues() operation. 1104 1105 You MUST have set the correct numbers of nonzeros per row in the call to 1106 MatCreateSeqSBAIJ(). 1107 1108 MUST be called before any calls to MatSetValues() 1109 1110 .seealso: MatCreateSeqSBAIJ 1111 @*/ 1112 int MatSeqSBAIJSetColumnIndices(Mat mat,int *indices) 1113 { 1114 int ierr,(*f)(Mat,int *); 1115 1116 PetscFunctionBegin; 1117 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1118 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void **)&f);CHKERRQ(ierr); 1119 if (f) { 1120 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1121 } else { 1122 SETERRQ(1,"Wrong type of matrix to set column indices"); 1123 } 1124 PetscFunctionReturn(0); 1125 } 1126 1127 #undef __FUNC__ 1128 #define __FUNC__ /*<a name="MatSetUpPreallocation_SeqSBAIJ"></a>*/"MatSetUpPreallocation_SeqSBAIJ" 1129 int MatSetUpPreallocation_SeqSBAIJ(Mat A) 1130 { 1131 int ierr; 1132 1133 PetscFunctionBegin; 1134 ierr = MatSeqSBAIJSetPreallocation(A,1,PETSC_DEFAULT,0);CHKERRQ(ierr); 1135 PetscFunctionReturn(0); 1136 } 1137 1138 /* -------------------------------------------------------------------*/ 1139 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1140 MatGetRow_SeqSBAIJ, 1141 MatRestoreRow_SeqSBAIJ, 1142 MatMult_SeqSBAIJ_N, 1143 MatMultAdd_SeqSBAIJ_N, 1144 MatMultTranspose_SeqSBAIJ, 1145 MatMultTransposeAdd_SeqSBAIJ, 1146 MatSolve_SeqSBAIJ_N, 1147 0, 1148 0, 1149 0, 1150 0, 1151 MatCholeskyFactor_SeqSBAIJ, 1152 0, 1153 MatTranspose_SeqSBAIJ, 1154 MatGetInfo_SeqSBAIJ, 1155 MatEqual_SeqSBAIJ, 1156 MatGetDiagonal_SeqSBAIJ, 1157 MatDiagonalScale_SeqSBAIJ, 1158 MatNorm_SeqSBAIJ, 1159 0, 1160 MatAssemblyEnd_SeqSBAIJ, 1161 0, 1162 MatSetOption_SeqSBAIJ, 1163 MatZeroEntries_SeqSBAIJ, 1164 MatZeroRows_SeqSBAIJ, 1165 0, 1166 0, 1167 MatCholeskyFactorSymbolic_SeqSBAIJ, 1168 MatCholeskyFactorNumeric_SeqSBAIJ_N, 1169 MatSetUpPreallocation_SeqSBAIJ, 1170 0, 1171 MatGetOwnershipRange_SeqSBAIJ, 1172 0, 1173 MatIncompleteCholeskyFactorSymbolic_SeqSBAIJ, 1174 0, 1175 0, 1176 MatDuplicate_SeqSBAIJ, 1177 0, 1178 0, 1179 0, 1180 0, 1181 0, 1182 MatGetSubMatrices_SeqSBAIJ, 1183 MatIncreaseOverlap_SeqSBAIJ, 1184 MatGetValues_SeqSBAIJ, 1185 0, 1186 MatPrintHelp_SeqSBAIJ, 1187 MatScale_SeqSBAIJ, 1188 0, 1189 0, 1190 0, 1191 MatGetBlockSize_SeqSBAIJ, 1192 MatGetRowIJ_SeqSBAIJ, 1193 MatRestoreRowIJ_SeqSBAIJ, 1194 0, 1195 0, 1196 0, 1197 0, 1198 0, 1199 0, 1200 MatSetValuesBlocked_SeqSBAIJ, 1201 MatGetSubMatrix_SeqSBAIJ, 1202 0, 1203 0, 1204 MatGetMaps_Petsc, 1205 0, 1206 0, 1207 0, 1208 0, 1209 0, 1210 0, 1211 MatGetRowMax_SeqSBAIJ}; 1212 1213 EXTERN_C_BEGIN 1214 #undef __FUNC__ 1215 #define __FUNC__ "MatStoreValues_SeqSBAIJ" 1216 int MatStoreValues_SeqSBAIJ(Mat mat) 1217 { 1218 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1219 int nz = aij->i[mat->m]*aij->bs*aij->bs2; 1220 int ierr; 1221 1222 PetscFunctionBegin; 1223 if (aij->nonew != 1) { 1224 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1225 } 1226 1227 /* allocate space for values if not already there */ 1228 if (!aij->saved_values) { 1229 aij->saved_values = (Scalar*)PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 1230 } 1231 1232 /* copy values over */ 1233 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 1234 PetscFunctionReturn(0); 1235 } 1236 EXTERN_C_END 1237 1238 EXTERN_C_BEGIN 1239 #undef __FUNC__ 1240 #define __FUNC__ "MatRetrieveValues_SeqSBAIJ" 1241 int MatRetrieveValues_SeqSBAIJ(Mat mat) 1242 { 1243 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1244 int nz = aij->i[mat->m]*aij->bs*aij->bs2,ierr; 1245 1246 PetscFunctionBegin; 1247 if (aij->nonew != 1) { 1248 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1249 } 1250 if (!aij->saved_values) { 1251 SETERRQ(1,"Must call MatStoreValues(A);first"); 1252 } 1253 1254 /* copy values over */ 1255 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 1256 PetscFunctionReturn(0); 1257 } 1258 EXTERN_C_END 1259 1260 EXTERN_C_BEGIN 1261 #undef __FUNC__ 1262 #define __FUNC__ "MatCreate_SeqSBAIJ" 1263 int MatCreate_SeqSBAIJ(Mat B) 1264 { 1265 Mat_SeqSBAIJ *b; 1266 int ierr,size; 1267 1268 PetscFunctionBegin; 1269 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1270 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1271 B->m = B->M = PetscMax(B->m,B->M); 1272 B->n = B->N = PetscMax(B->n,B->N); 1273 1274 B->data = (void*)(b = PetscNew(Mat_SeqSBAIJ));CHKPTRQ(b); 1275 ierr = PetscMemzero(b,sizeof(Mat_SeqSBAIJ));CHKERRQ(ierr); 1276 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1277 B->ops->destroy = MatDestroy_SeqSBAIJ; 1278 B->ops->view = MatView_SeqSBAIJ; 1279 B->factor = 0; 1280 B->lupivotthreshold = 1.0; 1281 B->mapping = 0; 1282 b->row = 0; 1283 b->icol = 0; 1284 b->reallocs = 0; 1285 b->saved_values = 0; 1286 1287 ierr = MapCreateMPI(B->comm,B->m,B->M,&B->rmap);CHKERRQ(ierr); 1288 ierr = MapCreateMPI(B->comm,B->n,B->N,&B->cmap);CHKERRQ(ierr); 1289 1290 b->sorted = PETSC_FALSE; 1291 b->roworiented = PETSC_TRUE; 1292 b->nonew = 0; 1293 b->diag = 0; 1294 b->solve_work = 0; 1295 b->mult_work = 0; 1296 b->spptr = 0; 1297 b->keepzeroedrows = PETSC_FALSE; 1298 1299 b->inew = 0; 1300 b->jnew = 0; 1301 b->anew = 0; 1302 b->a2anew = 0; 1303 b->permute = PETSC_FALSE; 1304 1305 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1306 "MatStoreValues_SeqSBAIJ", 1307 (void*)MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1308 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1309 "MatRetrieveValues_SeqSBAIJ", 1310 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1311 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1312 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1313 (void*)MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1314 PetscFunctionReturn(0); 1315 } 1316 EXTERN_C_END 1317 1318 #undef __FUNC__ 1319 #define __FUNC__ "MatSeqSBAIJSetPreallocation" 1320 /*@C 1321 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1322 compressed row) format. For good matrix assembly performance the 1323 user should preallocate the matrix storage by setting the parameter nz 1324 (or the array nnz). By setting these parameters accurately, performance 1325 during matrix assembly can be increased by more than a factor of 50. 1326 1327 Collective on Mat 1328 1329 Input Parameters: 1330 + A - the symmetric matrix 1331 . bs - size of block 1332 . nz - number of block nonzeros per block row (same for all rows) 1333 - nnz - array containing the number of block nonzeros in the various block rows 1334 (possibly different for each block row) or PETSC_NULL 1335 1336 Options Database Keys: 1337 . -mat_no_unroll - uses code that does not unroll the loops in the 1338 block calculations (much slower) 1339 . -mat_block_size - size of the blocks to use 1340 1341 Level: intermediate 1342 1343 Notes: 1344 The block AIJ format is compatible with standard Fortran 77 1345 storage. That is, the stored row and column indices can begin at 1346 either one (as in Fortran) or zero. See the users' manual for details. 1347 1348 Specify the preallocated storage with either nz or nnz (not both). 1349 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1350 allocation. For additional details, see the users manual chapter on 1351 matrices. 1352 1353 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1354 @*/ 1355 int MatSeqSBAIJSetPreallocation(Mat B,int bs,int nz,int *nnz) 1356 { 1357 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1358 int i,len,ierr,mbs,bs2; 1359 PetscTruth flg; 1360 int s_nz; 1361 1362 PetscFunctionBegin; 1363 B->preallocated = PETSC_TRUE; 1364 ierr = OptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1365 mbs = B->m/bs; 1366 bs2 = bs*bs; 1367 1368 if (mbs*bs != B->m) { 1369 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1370 } 1371 1372 if (nz < -2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than -2: value %d",nz); 1373 if (nnz) { 1374 for (i=0; i<mbs; i++) { 1375 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 1376 } 1377 } 1378 1379 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1380 if (!flg) { 1381 switch (bs) { 1382 case 1: 1383 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1384 B->ops->solve = MatSolve_SeqSBAIJ_1; 1385 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_1; 1386 B->ops->mult = MatMult_SeqSBAIJ_1; 1387 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1388 break; 1389 case 2: 1390 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1391 B->ops->solve = MatSolve_SeqSBAIJ_2; 1392 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_2; 1393 B->ops->mult = MatMult_SeqSBAIJ_2; 1394 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1395 break; 1396 case 3: 1397 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1398 B->ops->solve = MatSolve_SeqSBAIJ_3; 1399 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_3; 1400 B->ops->mult = MatMult_SeqSBAIJ_3; 1401 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1402 break; 1403 case 4: 1404 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1405 B->ops->solve = MatSolve_SeqSBAIJ_4; 1406 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_4; 1407 B->ops->mult = MatMult_SeqSBAIJ_4; 1408 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1409 break; 1410 case 5: 1411 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1412 B->ops->solve = MatSolve_SeqSBAIJ_5; 1413 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_5; 1414 B->ops->mult = MatMult_SeqSBAIJ_5; 1415 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1416 break; 1417 case 6: 1418 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1419 B->ops->solve = MatSolve_SeqSBAIJ_6; 1420 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_6; 1421 B->ops->mult = MatMult_SeqSBAIJ_6; 1422 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1423 break; 1424 case 7: 1425 B->ops->mult = MatMult_SeqSBAIJ_7; 1426 B->ops->solve = MatSolve_SeqSBAIJ_7; 1427 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_7; 1428 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1429 break; 1430 } 1431 } 1432 1433 b->mbs = mbs; 1434 b->nbs = mbs; 1435 b->imax = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(b->imax); 1436 if (!nnz) { 1437 if (nz == PETSC_DEFAULT) nz = 5; 1438 else if (nz <= 0) nz = 1; 1439 for (i=0; i<mbs; i++) { 1440 b->imax[i] = nz; 1441 } 1442 nz = nz*mbs; /* total nz */ 1443 } else { 1444 nz = 0; 1445 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1446 } 1447 /* s_nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1448 s_nz = nz; 1449 1450 /* allocate the matrix space */ 1451 len = s_nz*sizeof(int) + s_nz*bs2*sizeof(MatScalar) + (B->m+1)*sizeof(int); 1452 b->a = (MatScalar*)PetscMalloc(len);CHKPTRQ(b->a); 1453 ierr = PetscMemzero(b->a,s_nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1454 b->j = (int*)(b->a + s_nz*bs2); 1455 ierr = PetscMemzero(b->j,s_nz*sizeof(int));CHKERRQ(ierr); 1456 b->i = b->j + s_nz; 1457 b->singlemalloc = PETSC_TRUE; 1458 1459 /* pointer to beginning of each row */ 1460 b->i[0] = 0; 1461 for (i=1; i<mbs+1; i++) { 1462 b->i[i] = b->i[i-1] + b->imax[i-1]; 1463 } 1464 1465 /* b->ilen will count nonzeros in each block row so far. */ 1466 b->ilen = (int*)PetscMalloc((mbs+1)*sizeof(int)); 1467 PLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1468 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1469 1470 b->bs = bs; 1471 b->bs2 = bs2; 1472 b->s_nz = 0; 1473 b->s_maxnz = s_nz*bs2; 1474 1475 b->inew = 0; 1476 b->jnew = 0; 1477 b->anew = 0; 1478 b->a2anew = 0; 1479 b->permute = PETSC_FALSE; 1480 PetscFunctionReturn(0); 1481 } 1482 1483 1484 #undef __FUNC__ 1485 #define __FUNC__ "MatCreateSeqSBAIJ" 1486 /*@C 1487 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1488 compressed row) format. For good matrix assembly performance the 1489 user should preallocate the matrix storage by setting the parameter nz 1490 (or the array nnz). By setting these parameters accurately, performance 1491 during matrix assembly can be increased by more than a factor of 50. 1492 1493 Collective on MPI_Comm 1494 1495 Input Parameters: 1496 + comm - MPI communicator, set to PETSC_COMM_SELF 1497 . bs - size of block 1498 . m - number of rows, or number of columns 1499 . nz - number of block nonzeros per block row (same for all rows) 1500 - nnz - array containing the number of block nonzeros in the various block rows 1501 (possibly different for each block row) or PETSC_NULL 1502 1503 Output Parameter: 1504 . A - the symmetric matrix 1505 1506 Options Database Keys: 1507 . -mat_no_unroll - uses code that does not unroll the loops in the 1508 block calculations (much slower) 1509 . -mat_block_size - size of the blocks to use 1510 1511 Level: intermediate 1512 1513 Notes: 1514 The block AIJ format is compatible with standard Fortran 77 1515 storage. That is, the stored row and column indices can begin at 1516 either one (as in Fortran) or zero. See the users' manual for details. 1517 1518 Specify the preallocated storage with either nz or nnz (not both). 1519 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1520 allocation. For additional details, see the users manual chapter on 1521 matrices. 1522 1523 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1524 @*/ 1525 int MatCreateSeqSBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1526 { 1527 int ierr; 1528 1529 PetscFunctionBegin; 1530 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1531 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1532 ierr = MatSeqSBAIJSetPreallocation(*A,bs,nz,nnz);CHKERRQ(ierr); 1533 PetscFunctionReturn(0); 1534 } 1535 1536 #undef __FUNC__ 1537 #define __FUNC__ "MatDuplicate_SeqSBAIJ" 1538 int MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1539 { 1540 Mat C; 1541 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1542 int i,len,mbs = a->mbs,nz = a->s_nz,bs2 =a->bs2,ierr; 1543 1544 PetscFunctionBegin; 1545 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1546 1547 *B = 0; 1548 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 1549 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 1550 c = (Mat_SeqSBAIJ*)C->data; 1551 1552 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1553 C->preallocated = PETSC_TRUE; 1554 C->factor = A->factor; 1555 c->row = 0; 1556 c->icol = 0; 1557 c->saved_values = 0; 1558 c->keepzeroedrows = a->keepzeroedrows; 1559 C->assembled = PETSC_TRUE; 1560 1561 c->bs = a->bs; 1562 c->bs2 = a->bs2; 1563 c->mbs = a->mbs; 1564 c->nbs = a->nbs; 1565 1566 c->imax = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->imax); 1567 c->ilen = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->ilen); 1568 for (i=0; i<mbs; i++) { 1569 c->imax[i] = a->imax[i]; 1570 c->ilen[i] = a->ilen[i]; 1571 } 1572 1573 /* allocate the matrix space */ 1574 c->singlemalloc = PETSC_TRUE; 1575 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1576 c->a = (MatScalar*)PetscMalloc(len);CHKPTRQ(c->a); 1577 c->j = (int*)(c->a + nz*bs2); 1578 c->i = c->j + nz; 1579 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));CHKERRQ(ierr); 1580 if (mbs > 0) { 1581 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(int));CHKERRQ(ierr); 1582 if (cpvalues == MAT_COPY_VALUES) { 1583 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1584 } else { 1585 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1586 } 1587 } 1588 1589 PLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1590 c->sorted = a->sorted; 1591 c->roworiented = a->roworiented; 1592 c->nonew = a->nonew; 1593 1594 if (a->diag) { 1595 c->diag = (int*)PetscMalloc((mbs+1)*sizeof(int));CHKPTRQ(c->diag); 1596 PLogObjectMemory(C,(mbs+1)*sizeof(int)); 1597 for (i=0; i<mbs; i++) { 1598 c->diag[i] = a->diag[i]; 1599 } 1600 } else c->diag = 0; 1601 c->s_nz = a->s_nz; 1602 c->s_maxnz = a->s_maxnz; 1603 c->solve_work = 0; 1604 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1605 c->mult_work = 0; 1606 *B = C; 1607 ierr = FListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1608 PetscFunctionReturn(0); 1609 } 1610 1611 EXTERN_C_BEGIN 1612 #undef __FUNC__ 1613 #define __FUNC__ "MatLoad_SeqSBAIJ" 1614 int MatLoad_SeqSBAIJ(Viewer viewer,MatType type,Mat *A) 1615 { 1616 Mat_SeqSBAIJ *a; 1617 Mat B; 1618 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1; 1619 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,*s_browlengths,maskcount; 1620 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1621 int *masked,nmask,tmp,bs2,ishift; 1622 Scalar *aa; 1623 MPI_Comm comm = ((PetscObject)viewer)->comm; 1624 1625 PetscFunctionBegin; 1626 ierr = OptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1627 bs2 = bs*bs; 1628 1629 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1630 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1631 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1632 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1633 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1634 M = header[1]; N = header[2]; nz = header[3]; 1635 1636 if (header[3] < 0) { 1637 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 1638 } 1639 1640 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1641 1642 /* 1643 This code adds extra rows to make sure the number of rows is 1644 divisible by the blocksize 1645 */ 1646 mbs = M/bs; 1647 extra_rows = bs - M + bs*(mbs); 1648 if (extra_rows == bs) extra_rows = 0; 1649 else mbs++; 1650 if (extra_rows) { 1651 PLogInfo(0,"MatLoad_SeqSBAIJ:Padding loaded matrix to match blocksize\n"); 1652 } 1653 1654 /* read in row lengths */ 1655 rowlengths = (int*)PetscMalloc((M+extra_rows)*sizeof(int));CHKPTRQ(rowlengths); 1656 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1657 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1658 1659 /* read in column indices */ 1660 jj = (int*)PetscMalloc((nz+extra_rows)*sizeof(int));CHKPTRQ(jj); 1661 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1662 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1663 1664 /* loop over row lengths determining block row lengths */ 1665 browlengths = (int*)PetscMalloc(mbs*sizeof(int));CHKPTRQ(browlengths); 1666 s_browlengths = (int*)PetscMalloc(mbs*sizeof(int));CHKPTRQ(s_browlengths); 1667 ierr = PetscMemzero(s_browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1668 mask = (int*)PetscMalloc(2*mbs*sizeof(int));CHKPTRQ(mask); 1669 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1670 masked = mask + mbs; 1671 rowcount = 0; nzcount = 0; 1672 for (i=0; i<mbs; i++) { 1673 nmask = 0; 1674 for (j=0; j<bs; j++) { 1675 kmax = rowlengths[rowcount]; 1676 for (k=0; k<kmax; k++) { 1677 tmp = jj[nzcount++]/bs; /* block col. index */ 1678 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 1679 } 1680 rowcount++; 1681 } 1682 s_browlengths[i] += nmask; 1683 browlengths[i] = 2*s_browlengths[i]; 1684 1685 /* zero out the mask elements we set */ 1686 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1687 } 1688 1689 /* create our matrix */ 1690 ierr = MatCreateSeqSBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1691 B = *A; 1692 a = (Mat_SeqSBAIJ*)B->data; 1693 1694 /* set matrix "i" values */ 1695 a->i[0] = 0; 1696 for (i=1; i<= mbs; i++) { 1697 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 1698 a->ilen[i-1] = s_browlengths[i-1]; 1699 } 1700 a->s_nz = 0; 1701 for (i=0; i<mbs; i++) a->s_nz += s_browlengths[i]; 1702 1703 /* read in nonzero values */ 1704 aa = (Scalar*)PetscMalloc((nz+extra_rows)*sizeof(Scalar));CHKPTRQ(aa); 1705 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1706 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1707 1708 /* set "a" and "j" values into matrix */ 1709 nzcount = 0; jcount = 0; 1710 for (i=0; i<mbs; i++) { 1711 nzcountb = nzcount; 1712 nmask = 0; 1713 for (j=0; j<bs; j++) { 1714 kmax = rowlengths[i*bs+j]; 1715 for (k=0; k<kmax; k++) { 1716 tmp = jj[nzcount++]/bs; /* block col. index */ 1717 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 1718 } 1719 rowcount++; 1720 } 1721 /* sort the masked values */ 1722 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1723 1724 /* set "j" values into matrix */ 1725 maskcount = 1; 1726 for (j=0; j<nmask; j++) { 1727 a->j[jcount++] = masked[j]; 1728 mask[masked[j]] = maskcount++; 1729 } 1730 1731 /* set "a" values into matrix */ 1732 ishift = bs2*a->i[i]; 1733 for (j=0; j<bs; j++) { 1734 kmax = rowlengths[i*bs+j]; 1735 for (k=0; k<kmax; k++) { 1736 tmp = jj[nzcountb]/bs ; /* block col. index */ 1737 if (tmp >= i){ 1738 block = mask[tmp] - 1; 1739 point = jj[nzcountb] - bs*tmp; 1740 idx = ishift + bs2*block + j + bs*point; 1741 a->a[idx] = aa[nzcountb]; 1742 } 1743 nzcountb++; 1744 } 1745 } 1746 /* zero out the mask elements we set */ 1747 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1748 } 1749 if (jcount != a->s_nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 1750 1751 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1752 ierr = PetscFree(browlengths);CHKERRQ(ierr); 1753 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 1754 ierr = PetscFree(aa);CHKERRQ(ierr); 1755 ierr = PetscFree(jj);CHKERRQ(ierr); 1756 ierr = PetscFree(mask);CHKERRQ(ierr); 1757 1758 B->assembled = PETSC_TRUE; 1759 ierr = MatView_Private(B);CHKERRQ(ierr); 1760 PetscFunctionReturn(0); 1761 } 1762 EXTERN_C_END 1763 1764 1765 1766