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