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