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