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