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