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