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