1 #define PETSCMAT_DLL 2 3 /* 4 Defines the basic matrix operations for the SBAIJ (compressed row) 5 matrix storage format. 6 */ 7 #include "src/mat/impls/baij/seq/baij.h" /*I "petscmat.h" I*/ 8 #include "src/inline/spops.h" 9 #include "src/mat/impls/sbaij/seq/sbaij.h" 10 11 #define CHUNKSIZE 10 12 13 /* 14 Checks for missing diagonals 15 */ 16 #undef __FUNCT__ 17 #define __FUNCT__ "MatMissingDiagonal_SeqSBAIJ" 18 PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A) 19 { 20 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 21 PetscErrorCode ierr; 22 PetscInt *diag,*jj = a->j,i; 23 24 PetscFunctionBegin; 25 ierr = MatMarkDiagonal_SeqSBAIJ(A);CHKERRQ(ierr); 26 diag = a->diag; 27 for (i=0; i<a->mbs; i++) { 28 if (jj[diag[i]] != i) SETERRQ1(PETSC_ERR_ARG_CORRUPT,"Matrix is missing diagonal number %D",i); 29 } 30 PetscFunctionReturn(0); 31 } 32 33 #undef __FUNCT__ 34 #define __FUNCT__ "MatMarkDiagonal_SeqSBAIJ" 35 PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A) 36 { 37 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 38 PetscErrorCode ierr; 39 PetscInt i; 40 41 PetscFunctionBegin; 42 if (!a->diag) { 43 ierr = PetscMalloc(a->mbs*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 44 } 45 for (i=0; i<a->mbs; i++) a->diag[i] = a->i[i]; 46 PetscFunctionReturn(0); 47 } 48 49 #undef __FUNCT__ 50 #define __FUNCT__ "MatGetRowIJ_SeqSBAIJ" 51 static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 52 { 53 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 54 PetscInt n = a->mbs,i; 55 56 PetscFunctionBegin; 57 *nn = n; 58 if (!ia) PetscFunctionReturn(0); 59 60 if (oshift == 1) { 61 /* temporarily add 1 to i and j indices */ 62 PetscInt nz = a->i[n]; 63 for (i=0; i<nz; i++) a->j[i]++; 64 for (i=0; i<n+1; i++) a->i[i]++; 65 *ia = a->i; *ja = a->j; 66 } else { 67 *ia = a->i; *ja = a->j; 68 } 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "MatRestoreRowIJ_SeqSBAIJ" 74 static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 75 { 76 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 77 PetscInt i,n = a->mbs; 78 79 PetscFunctionBegin; 80 if (!ia) PetscFunctionReturn(0); 81 82 if (oshift == 1) { 83 PetscInt nz = a->i[n]-1; 84 for (i=0; i<nz; i++) a->j[i]--; 85 for (i=0; i<n+1; i++) a->i[i]--; 86 } 87 PetscFunctionReturn(0); 88 } 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "MatDestroy_SeqSBAIJ" 92 PetscErrorCode MatDestroy_SeqSBAIJ(Mat A) 93 { 94 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 95 PetscErrorCode ierr; 96 97 PetscFunctionBegin; 98 #if defined(PETSC_USE_LOG) 99 PetscLogObjectState((PetscObject)A,"Rows=%D, NZ=%D",A->rmap.N,a->nz); 100 #endif 101 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 102 if (a->row) { 103 ierr = ISDestroy(a->row);CHKERRQ(ierr); 104 } 105 ierr = PetscFree(a->diag);CHKERRQ(ierr); 106 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 107 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 108 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 109 ierr = PetscFree(a->relax_work);CHKERRQ(ierr); 110 ierr = PetscFree(a->solves_work);CHKERRQ(ierr); 111 ierr = PetscFree(a->mult_work);CHKERRQ(ierr); 112 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 113 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 114 115 ierr = PetscFree(a->inew);CHKERRQ(ierr); 116 ierr = PetscFree(a);CHKERRQ(ierr); 117 118 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 119 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 120 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 121 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C","",PETSC_NULL);CHKERRQ(ierr); 122 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 123 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 127 #undef __FUNCT__ 128 #define __FUNCT__ "MatSetOption_SeqSBAIJ" 129 PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op) 130 { 131 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 switch (op) { 136 case MAT_ROW_ORIENTED: 137 a->roworiented = PETSC_TRUE; 138 break; 139 case MAT_COLUMN_ORIENTED: 140 a->roworiented = PETSC_FALSE; 141 break; 142 case MAT_COLUMNS_SORTED: 143 a->sorted = PETSC_TRUE; 144 break; 145 case MAT_COLUMNS_UNSORTED: 146 a->sorted = PETSC_FALSE; 147 break; 148 case MAT_KEEP_ZEROED_ROWS: 149 a->keepzeroedrows = PETSC_TRUE; 150 break; 151 case MAT_NO_NEW_NONZERO_LOCATIONS: 152 a->nonew = 1; 153 break; 154 case MAT_NEW_NONZERO_LOCATION_ERR: 155 a->nonew = -1; 156 break; 157 case MAT_NEW_NONZERO_ALLOCATION_ERR: 158 a->nonew = -2; 159 break; 160 case MAT_YES_NEW_NONZERO_LOCATIONS: 161 a->nonew = 0; 162 break; 163 case MAT_ROWS_SORTED: 164 case MAT_ROWS_UNSORTED: 165 case MAT_YES_NEW_DIAGONALS: 166 case MAT_IGNORE_OFF_PROC_ENTRIES: 167 case MAT_USE_HASH_TABLE: 168 ierr = PetscInfo1(A,"Option %d ignored\n",op);CHKERRQ(ierr); 169 break; 170 case MAT_NO_NEW_DIAGONALS: 171 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 172 case MAT_NOT_SYMMETRIC: 173 case MAT_NOT_STRUCTURALLY_SYMMETRIC: 174 case MAT_HERMITIAN: 175 SETERRQ(PETSC_ERR_SUP,"Matrix must be symmetric"); 176 case MAT_SYMMETRIC: 177 case MAT_STRUCTURALLY_SYMMETRIC: 178 case MAT_NOT_HERMITIAN: 179 case MAT_SYMMETRY_ETERNAL: 180 case MAT_NOT_SYMMETRY_ETERNAL: 181 case MAT_IGNORE_LOWER_TRIANGULAR: 182 a->ignore_ltriangular = PETSC_TRUE; 183 break; 184 case MAT_ERROR_LOWER_TRIANGULAR: 185 a->ignore_ltriangular = PETSC_FALSE; 186 break; 187 case MAT_GETROW_UPPERTRIANGULAR: 188 a->getrow_utriangular = PETSC_TRUE; 189 break; 190 default: 191 SETERRQ1(PETSC_ERR_SUP,"unknown option %d",op); 192 } 193 PetscFunctionReturn(0); 194 } 195 196 #undef __FUNCT__ 197 #define __FUNCT__ "MatGetRow_SeqSBAIJ" 198 PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **v) 199 { 200 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 201 PetscErrorCode ierr; 202 PetscInt itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*cols_i,bs2; 203 MatScalar *aa,*aa_i; 204 PetscScalar *v_i; 205 206 PetscFunctionBegin; 207 if (A && !a->getrow_utriangular) SETERRQ(PETSC_ERR_SUP,"MatGetRow is not supported for SBAIJ matrix format. Getting the upper triangular part of row, run with -mat_getrow_uppertriangular, call MatSetOption(mat,MAT_GETROW_UPPERTRIANGULAR) or MatGetRowUpperTriangular()"); 208 /* Get the upper triangular part of the row */ 209 bs = A->rmap.bs; 210 ai = a->i; 211 aj = a->j; 212 aa = a->a; 213 bs2 = a->bs2; 214 215 if (row < 0 || row >= A->rmap.N) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Row %D out of range", row); 216 217 bn = row/bs; /* Block number */ 218 bp = row % bs; /* Block position */ 219 M = ai[bn+1] - ai[bn]; 220 *ncols = bs*M; 221 222 if (v) { 223 *v = 0; 224 if (*ncols) { 225 ierr = PetscMalloc((*ncols+row)*sizeof(PetscScalar),v);CHKERRQ(ierr); 226 for (i=0; i<M; i++) { /* for each block in the block row */ 227 v_i = *v + i*bs; 228 aa_i = aa + bs2*(ai[bn] + i); 229 for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];} 230 } 231 } 232 } 233 234 if (cols) { 235 *cols = 0; 236 if (*ncols) { 237 ierr = PetscMalloc((*ncols+row)*sizeof(PetscInt),cols);CHKERRQ(ierr); 238 for (i=0; i<M; i++) { /* for each block in the block row */ 239 cols_i = *cols + i*bs; 240 itmp = bs*aj[ai[bn] + i]; 241 for (j=0; j<bs; j++) {cols_i[j] = itmp++;} 242 } 243 } 244 } 245 246 /*search column A(0:row-1,row) (=A(row,0:row-1)). Could be expensive! */ 247 /* this segment is currently removed, so only entries in the upper triangle are obtained */ 248 #ifdef column_search 249 v_i = *v + M*bs; 250 cols_i = *cols + M*bs; 251 for (i=0; i<bn; i++){ /* for each block row */ 252 M = ai[i+1] - ai[i]; 253 for (j=0; j<M; j++){ 254 itmp = aj[ai[i] + j]; /* block column value */ 255 if (itmp == bn){ 256 aa_i = aa + bs2*(ai[i] + j) + bs*bp; 257 for (k=0; k<bs; k++) { 258 *cols_i++ = i*bs+k; 259 *v_i++ = aa_i[k]; 260 } 261 *ncols += bs; 262 break; 263 } 264 } 265 } 266 #endif 267 PetscFunctionReturn(0); 268 } 269 270 #undef __FUNCT__ 271 #define __FUNCT__ "MatRestoreRow_SeqSBAIJ" 272 PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 273 { 274 PetscErrorCode ierr; 275 276 PetscFunctionBegin; 277 if (idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);} 278 if (v) {ierr = PetscFree(*v);CHKERRQ(ierr);} 279 PetscFunctionReturn(0); 280 } 281 282 #undef __FUNCT__ 283 #define __FUNCT__ "MatGetRowUpperTriangular_SeqSBAIJ" 284 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A) 285 { 286 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 287 288 PetscFunctionBegin; 289 a->getrow_utriangular = PETSC_TRUE; 290 PetscFunctionReturn(0); 291 } 292 #undef __FUNCT__ 293 #define __FUNCT__ "MatRestoreRowUpperTriangular_SeqSBAIJ" 294 PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A) 295 { 296 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 297 298 PetscFunctionBegin; 299 a->getrow_utriangular = PETSC_FALSE; 300 PetscFunctionReturn(0); 301 } 302 303 #undef __FUNCT__ 304 #define __FUNCT__ "MatTranspose_SeqSBAIJ" 305 PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,Mat *B) 306 { 307 PetscErrorCode ierr; 308 PetscFunctionBegin; 309 ierr = MatDuplicate(A,MAT_COPY_VALUES,B);CHKERRQ(ierr); 310 PetscFunctionReturn(0); 311 } 312 313 #undef __FUNCT__ 314 #define __FUNCT__ "MatView_SeqSBAIJ_ASCII" 315 static PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 316 { 317 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 318 PetscErrorCode ierr; 319 PetscInt i,j,bs = A->rmap.bs,k,l,bs2=a->bs2; 320 const char *name; 321 PetscViewerFormat format; 322 323 PetscFunctionBegin; 324 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 325 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 326 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 327 ierr = PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);CHKERRQ(ierr); 328 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 329 SETERRQ(PETSC_ERR_SUP,"Matlab format not supported"); 330 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 331 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 332 for (i=0; i<a->mbs; i++) { 333 for (j=0; j<bs; j++) { 334 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 335 for (k=a->i[i]; k<a->i[i+1]; k++) { 336 for (l=0; l<bs; l++) { 337 #if defined(PETSC_USE_COMPLEX) 338 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 339 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 340 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 341 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 342 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 343 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 344 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 345 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 346 } 347 #else 348 if (a->a[bs2*k + l*bs + j] != 0.0) { 349 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 350 } 351 #endif 352 } 353 } 354 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 355 } 356 } 357 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 358 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 359 PetscFunctionReturn(0); 360 } else { 361 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 362 for (i=0; i<a->mbs; i++) { 363 for (j=0; j<bs; j++) { 364 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 365 for (k=a->i[i]; k<a->i[i+1]; k++) { 366 for (l=0; l<bs; l++) { 367 #if defined(PETSC_USE_COMPLEX) 368 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 369 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 370 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 371 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 372 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 373 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 374 } else { 375 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 376 } 377 #else 378 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 379 #endif 380 } 381 } 382 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 383 } 384 } 385 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 386 } 387 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 388 PetscFunctionReturn(0); 389 } 390 391 #undef __FUNCT__ 392 #define __FUNCT__ "MatView_SeqSBAIJ_Draw_Zoom" 393 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 394 { 395 Mat A = (Mat) Aa; 396 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 397 PetscErrorCode ierr; 398 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap.bs,bs2=a->bs2; 399 PetscMPIInt rank; 400 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 401 MatScalar *aa; 402 MPI_Comm comm; 403 PetscViewer viewer; 404 405 PetscFunctionBegin; 406 /* 407 This is nasty. If this is called from an originally parallel matrix 408 then all processes call this,but only the first has the matrix so the 409 rest should return immediately. 410 */ 411 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 412 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 413 if (rank) PetscFunctionReturn(0); 414 415 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 416 417 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 418 PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric"); 419 420 /* loop over matrix elements drawing boxes */ 421 color = PETSC_DRAW_BLUE; 422 for (i=0,row=0; i<mbs; i++,row+=bs) { 423 for (j=a->i[i]; j<a->i[i+1]; j++) { 424 y_l = A->rmap.N - row - 1.0; y_r = y_l + 1.0; 425 x_l = a->j[j]*bs; x_r = x_l + 1.0; 426 aa = a->a + j*bs2; 427 for (k=0; k<bs; k++) { 428 for (l=0; l<bs; l++) { 429 if (PetscRealPart(*aa++) >= 0.) continue; 430 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 431 } 432 } 433 } 434 } 435 color = PETSC_DRAW_CYAN; 436 for (i=0,row=0; i<mbs; i++,row+=bs) { 437 for (j=a->i[i]; j<a->i[i+1]; j++) { 438 y_l = A->rmap.N - row - 1.0; y_r = y_l + 1.0; 439 x_l = a->j[j]*bs; x_r = x_l + 1.0; 440 aa = a->a + j*bs2; 441 for (k=0; k<bs; k++) { 442 for (l=0; l<bs; l++) { 443 if (PetscRealPart(*aa++) != 0.) continue; 444 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 445 } 446 } 447 } 448 } 449 450 color = PETSC_DRAW_RED; 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->rmap.N - 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 = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 460 } 461 } 462 } 463 } 464 PetscFunctionReturn(0); 465 } 466 467 #undef __FUNCT__ 468 #define __FUNCT__ "MatView_SeqSBAIJ_Draw" 469 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 470 { 471 PetscErrorCode ierr; 472 PetscReal xl,yl,xr,yr,w,h; 473 PetscDraw draw; 474 PetscTruth isnull; 475 476 PetscFunctionBegin; 477 478 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 479 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 480 481 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 482 xr = A->rmap.N; yr = A->rmap.N; h = yr/10.0; w = xr/10.0; 483 xr += w; yr += h; xl = -w; yl = -h; 484 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 485 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 486 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 487 PetscFunctionReturn(0); 488 } 489 490 #undef __FUNCT__ 491 #define __FUNCT__ "MatView_SeqSBAIJ" 492 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 493 { 494 PetscErrorCode ierr; 495 PetscTruth iascii,isdraw; 496 497 PetscFunctionBegin; 498 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 499 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 500 if (iascii){ 501 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 502 } else if (isdraw) { 503 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 504 } else { 505 Mat B; 506 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 507 ierr = MatView(B,viewer);CHKERRQ(ierr); 508 ierr = MatDestroy(B);CHKERRQ(ierr); 509 } 510 PetscFunctionReturn(0); 511 } 512 513 514 #undef __FUNCT__ 515 #define __FUNCT__ "MatGetValues_SeqSBAIJ" 516 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 517 { 518 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 519 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 520 PetscInt *ai = a->i,*ailen = a->ilen; 521 PetscInt brow,bcol,ridx,cidx,bs=A->rmap.bs,bs2=a->bs2; 522 MatScalar *ap,*aa = a->a,zero = 0.0; 523 524 PetscFunctionBegin; 525 for (k=0; k<m; k++) { /* loop over rows */ 526 row = im[k]; brow = row/bs; 527 if (row < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); 528 if (row >= A->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap.N-1); 529 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 530 nrow = ailen[brow]; 531 for (l=0; l<n; l++) { /* loop over columns */ 532 if (in[l] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); 533 if (in[l] >= A->cmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap.n-1); 534 col = in[l] ; 535 bcol = col/bs; 536 cidx = col%bs; 537 ridx = row%bs; 538 high = nrow; 539 low = 0; /* assume unsorted */ 540 while (high-low > 5) { 541 t = (low+high)/2; 542 if (rp[t] > bcol) high = t; 543 else low = t; 544 } 545 for (i=low; i<high; i++) { 546 if (rp[i] > bcol) break; 547 if (rp[i] == bcol) { 548 *v++ = ap[bs2*i+bs*cidx+ridx]; 549 goto finished; 550 } 551 } 552 *v++ = zero; 553 finished:; 554 } 555 } 556 PetscFunctionReturn(0); 557 } 558 559 560 #undef __FUNCT__ 561 #define __FUNCT__ "MatSetValuesBlocked_SeqSBAIJ" 562 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 563 { 564 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 565 PetscErrorCode ierr; 566 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 567 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen; 568 PetscInt *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap.bs,stepval; 569 PetscTruth roworiented=a->roworiented; 570 const MatScalar *value = v; 571 MatScalar *ap,*aa = a->a,*bap; 572 573 PetscFunctionBegin; 574 if (roworiented) { 575 stepval = (n-1)*bs; 576 } else { 577 stepval = (m-1)*bs; 578 } 579 for (k=0; k<m; k++) { /* loop over added rows */ 580 row = im[k]; 581 if (row < 0) continue; 582 #if defined(PETSC_USE_DEBUG) 583 if (row >= a->mbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,a->mbs-1); 584 #endif 585 rp = aj + ai[row]; 586 ap = aa + bs2*ai[row]; 587 rmax = imax[row]; 588 nrow = ailen[row]; 589 low = 0; 590 high = nrow; 591 for (l=0; l<n; l++) { /* loop over added columns */ 592 if (in[l] < 0) continue; 593 col = in[l]; 594 #if defined(PETSC_USE_DEBUG) 595 if (col >= a->nbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",col,a->nbs-1); 596 #endif 597 if (col < row) continue; /* ignore lower triangular block */ 598 if (roworiented) { 599 value = v + k*(stepval+bs)*bs + l*bs; 600 } else { 601 value = v + l*(stepval+bs)*bs + k*bs; 602 } 603 if (col <= lastcol) low = 0; else high = nrow; 604 lastcol = col; 605 while (high-low > 7) { 606 t = (low+high)/2; 607 if (rp[t] > col) high = t; 608 else low = t; 609 } 610 for (i=low; i<high; i++) { 611 if (rp[i] > col) break; 612 if (rp[i] == col) { 613 bap = ap + bs2*i; 614 if (roworiented) { 615 if (is == ADD_VALUES) { 616 for (ii=0; ii<bs; ii++,value+=stepval) { 617 for (jj=ii; jj<bs2; jj+=bs) { 618 bap[jj] += *value++; 619 } 620 } 621 } else { 622 for (ii=0; ii<bs; ii++,value+=stepval) { 623 for (jj=ii; jj<bs2; jj+=bs) { 624 bap[jj] = *value++; 625 } 626 } 627 } 628 } else { 629 if (is == ADD_VALUES) { 630 for (ii=0; ii<bs; ii++,value+=stepval) { 631 for (jj=0; jj<bs; jj++) { 632 *bap++ += *value++; 633 } 634 } 635 } else { 636 for (ii=0; ii<bs; ii++,value+=stepval) { 637 for (jj=0; jj<bs; jj++) { 638 *bap++ = *value++; 639 } 640 } 641 } 642 } 643 goto noinsert2; 644 } 645 } 646 if (nonew == 1) goto noinsert2; 647 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 648 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew); 649 N = nrow++ - 1; high++; 650 /* shift up all the later entries in this row */ 651 for (ii=N; ii>=i; ii--) { 652 rp[ii+1] = rp[ii]; 653 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 654 } 655 if (N >= i) { 656 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 657 } 658 rp[i] = col; 659 bap = ap + bs2*i; 660 if (roworiented) { 661 for (ii=0; ii<bs; ii++,value+=stepval) { 662 for (jj=ii; jj<bs2; jj+=bs) { 663 bap[jj] = *value++; 664 } 665 } 666 } else { 667 for (ii=0; ii<bs; ii++,value+=stepval) { 668 for (jj=0; jj<bs; jj++) { 669 *bap++ = *value++; 670 } 671 } 672 } 673 noinsert2:; 674 low = i; 675 } 676 ailen[row] = nrow; 677 } 678 PetscFunctionReturn(0); 679 } 680 681 #undef __FUNCT__ 682 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 683 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 684 { 685 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 686 PetscErrorCode ierr; 687 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 688 PetscInt m = A->rmap.N,*ip,N,*ailen = a->ilen; 689 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 690 MatScalar *aa = a->a,*ap; 691 692 PetscFunctionBegin; 693 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 694 695 if (m) rmax = ailen[0]; 696 for (i=1; i<mbs; i++) { 697 /* move each row back by the amount of empty slots (fshift) before it*/ 698 fshift += imax[i-1] - ailen[i-1]; 699 rmax = PetscMax(rmax,ailen[i]); 700 if (fshift) { 701 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 702 N = ailen[i]; 703 for (j=0; j<N; j++) { 704 ip[j-fshift] = ip[j]; 705 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 706 } 707 } 708 ai[i] = ai[i-1] + ailen[i-1]; 709 } 710 if (mbs) { 711 fshift += imax[mbs-1] - ailen[mbs-1]; 712 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 713 } 714 /* reset ilen and imax for each row */ 715 for (i=0; i<mbs; i++) { 716 ailen[i] = imax[i] = ai[i+1] - ai[i]; 717 } 718 a->nz = ai[mbs]; 719 720 /* diagonals may have moved, reset it */ 721 if (a->diag) { 722 ierr = PetscMemcpy(a->diag,ai,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 723 } 724 ierr = PetscInfo5(A,"Matrix size: %D X %D, block size %D; storage space: %D unneeded, %D used\n",m,A->rmap.N,A->rmap.bs,fshift*bs2,a->nz*bs2);CHKERRQ(ierr); 725 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 726 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 727 a->reallocs = 0; 728 A->info.nz_unneeded = (PetscReal)fshift*bs2; 729 PetscFunctionReturn(0); 730 } 731 732 /* 733 This function returns an array of flags which indicate the locations of contiguous 734 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 735 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 736 Assume: sizes should be long enough to hold all the values. 737 */ 738 #undef __FUNCT__ 739 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 740 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 741 { 742 PetscInt i,j,k,row; 743 PetscTruth flg; 744 745 PetscFunctionBegin; 746 for (i=0,j=0; i<n; j++) { 747 row = idx[i]; 748 if (row%bs!=0) { /* Not the begining of a block */ 749 sizes[j] = 1; 750 i++; 751 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 752 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 753 i++; 754 } else { /* Begining of the block, so check if the complete block exists */ 755 flg = PETSC_TRUE; 756 for (k=1; k<bs; k++) { 757 if (row+k != idx[i+k]) { /* break in the block */ 758 flg = PETSC_FALSE; 759 break; 760 } 761 } 762 if (flg) { /* No break in the bs */ 763 sizes[j] = bs; 764 i+= bs; 765 } else { 766 sizes[j] = 1; 767 i++; 768 } 769 } 770 } 771 *bs_max = j; 772 PetscFunctionReturn(0); 773 } 774 775 776 /* Only add/insert a(i,j) with i<=j (blocks). 777 Any a(i,j) with i>j input by user is ingored. 778 */ 779 780 #undef __FUNCT__ 781 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 782 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 783 { 784 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 785 PetscErrorCode ierr; 786 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 787 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 788 PetscInt *aj=a->j,nonew=a->nonew,bs=A->rmap.bs,brow,bcol; 789 PetscInt ridx,cidx,bs2=a->bs2; 790 MatScalar *ap,value,*aa=a->a,*bap; 791 792 PetscFunctionBegin; 793 for (k=0; k<m; k++) { /* loop over added rows */ 794 row = im[k]; /* row number */ 795 brow = row/bs; /* block row number */ 796 if (row < 0) continue; 797 #if defined(PETSC_USE_DEBUG) 798 if (row >= A->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap.N-1); 799 #endif 800 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 801 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 802 rmax = imax[brow]; /* maximum space allocated for this row */ 803 nrow = ailen[brow]; /* actual length of this row */ 804 low = 0; 805 806 for (l=0; l<n; l++) { /* loop over added columns */ 807 if (in[l] < 0) continue; 808 #if defined(PETSC_USE_DEBUG) 809 if (in[l] >= A->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->rmap.N-1); 810 #endif 811 col = in[l]; 812 bcol = col/bs; /* block col number */ 813 814 if (brow > bcol) { 815 if (a->ignore_ltriangular){ 816 continue; /* ignore lower triangular values */ 817 } else { 818 SETERRQ(PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR)"); 819 } 820 } 821 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 (col <= lastcol) low = 0; high = nrow; 833 lastcol = col; 834 while (high-low > 7) { 835 t = (low+high)/2; 836 if (rp[t] > bcol) high = t; 837 else low = t; 838 } 839 for (i=low; i<high; 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 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 846 if (brow == bcol && ridx < cidx){ 847 bap = ap + bs2*i + bs*ridx + cidx; 848 if (is == ADD_VALUES) *bap += value; 849 else *bap = value; 850 } 851 goto noinsert1; 852 } 853 } 854 855 if (nonew == 1) goto noinsert1; 856 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 857 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew); 858 859 N = nrow++ - 1; high++; 860 /* shift up all the later entries in this row */ 861 for (ii=N; ii>=i; ii--) { 862 rp[ii+1] = rp[ii]; 863 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 864 } 865 if (N>=i) { 866 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 867 } 868 rp[i] = bcol; 869 ap[bs2*i + bs*cidx + ridx] = value; 870 noinsert1:; 871 low = i; 872 } 873 } /* end of loop over added columns */ 874 ailen[brow] = nrow; 875 } /* end of loop over added rows */ 876 PetscFunctionReturn(0); 877 } 878 879 EXTERN PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ(Mat,IS,MatFactorInfo*,Mat*); 880 EXTERN PetscErrorCode MatCholeskyFactor_SeqSBAIJ(Mat,IS,MatFactorInfo*); 881 EXTERN PetscErrorCode MatIncreaseOverlap_SeqSBAIJ(Mat,PetscInt,IS[],PetscInt); 882 EXTERN PetscErrorCode MatGetSubMatrix_SeqSBAIJ(Mat,IS,IS,PetscInt,MatReuse,Mat*); 883 EXTERN PetscErrorCode MatGetSubMatrices_SeqSBAIJ(Mat,PetscInt,const IS[],const IS[],MatReuse,Mat*[]); 884 EXTERN PetscErrorCode MatScale_SeqSBAIJ(Mat,PetscScalar); 885 EXTERN PetscErrorCode MatNorm_SeqSBAIJ(Mat,NormType,PetscReal *); 886 EXTERN PetscErrorCode MatEqual_SeqSBAIJ(Mat,Mat,PetscTruth*); 887 EXTERN PetscErrorCode MatGetDiagonal_SeqSBAIJ(Mat,Vec); 888 EXTERN PetscErrorCode MatDiagonalScale_SeqSBAIJ(Mat,Vec,Vec); 889 EXTERN PetscErrorCode MatGetInfo_SeqSBAIJ(Mat,MatInfoType,MatInfo *); 890 EXTERN PetscErrorCode MatZeroEntries_SeqSBAIJ(Mat); 891 EXTERN PetscErrorCode MatGetRowMax_SeqSBAIJ(Mat,Vec); 892 893 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_N(Mat,Vec,Vec); 894 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_1(Mat,Vec,Vec); 895 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_2(Mat,Vec,Vec); 896 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_3(Mat,Vec,Vec); 897 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_4(Mat,Vec,Vec); 898 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_5(Mat,Vec,Vec); 899 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_6(Mat,Vec,Vec); 900 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_7(Mat,Vec,Vec); 901 902 EXTERN PetscErrorCode MatSolves_SeqSBAIJ_1(Mat,Vecs,Vecs); 903 904 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_1_NaturalOrdering(Mat,Vec,Vec); 905 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_2_NaturalOrdering(Mat,Vec,Vec); 906 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_3_NaturalOrdering(Mat,Vec,Vec); 907 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_4_NaturalOrdering(Mat,Vec,Vec); 908 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_5_NaturalOrdering(Mat,Vec,Vec); 909 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_6_NaturalOrdering(Mat,Vec,Vec); 910 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_7_NaturalOrdering(Mat,Vec,Vec); 911 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_N_NaturalOrdering(Mat,Vec,Vec); 912 913 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat,MatFactorInfo*,Mat*); 914 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1(Mat,MatFactorInfo*,Mat*); 915 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat,MatFactorInfo*,Mat*); 916 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_3(Mat,MatFactorInfo*,Mat*); 917 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_4(Mat,MatFactorInfo*,Mat*); 918 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_5(Mat,MatFactorInfo*,Mat*); 919 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_6(Mat,MatFactorInfo*,Mat*); 920 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_7(Mat,MatFactorInfo*,Mat*); 921 EXTERN PetscErrorCode MatGetInertia_SeqSBAIJ(Mat,PetscInt*,PetscInt*,PetscInt*); 922 923 EXTERN PetscErrorCode MatMult_SeqSBAIJ_1(Mat,Vec,Vec); 924 EXTERN PetscErrorCode MatMult_SeqSBAIJ_2(Mat,Vec,Vec); 925 EXTERN PetscErrorCode MatMult_SeqSBAIJ_3(Mat,Vec,Vec); 926 EXTERN PetscErrorCode MatMult_SeqSBAIJ_4(Mat,Vec,Vec); 927 EXTERN PetscErrorCode MatMult_SeqSBAIJ_5(Mat,Vec,Vec); 928 EXTERN PetscErrorCode MatMult_SeqSBAIJ_6(Mat,Vec,Vec); 929 EXTERN PetscErrorCode MatMult_SeqSBAIJ_7(Mat,Vec,Vec); 930 EXTERN PetscErrorCode MatMult_SeqSBAIJ_N(Mat,Vec,Vec); 931 932 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_1(Mat,Vec,Vec,Vec); 933 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_2(Mat,Vec,Vec,Vec); 934 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_3(Mat,Vec,Vec,Vec); 935 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_4(Mat,Vec,Vec,Vec); 936 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_5(Mat,Vec,Vec,Vec); 937 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_6(Mat,Vec,Vec,Vec); 938 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_7(Mat,Vec,Vec,Vec); 939 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_N(Mat,Vec,Vec,Vec); 940 941 #ifdef HAVE_MatICCFactor 942 /* modified from MatILUFactor_SeqSBAIJ, needs further work! */ 943 #undef __FUNCT__ 944 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 945 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,MatFactorInfo *info) 946 { 947 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 948 Mat outA; 949 PetscErrorCode ierr; 950 PetscTruth row_identity,col_identity; 951 952 PetscFunctionBegin; 953 outA = inA; 954 inA->factor = FACTOR_CHOLESKY; 955 956 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 957 /* 958 Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver 959 for ILU(0) factorization with natural ordering 960 */ 961 switch (a->rmap.bs) { 962 case 1: 963 inA->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering; 964 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering; 965 inA->ops->solves = MatSolves_SeqSBAIJ_1; 966 ierr = PetscInfo((inA,"Using special in-place natural ordering solvetrans BS=1\n");CHKERRQ(ierr); 967 case 2: 968 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 969 inA->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering; 970 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering; 971 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=2\n");CHKERRQ(ierr); 972 break; 973 case 3: 974 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 975 inA->ops->solve = MatSolve_SeqSBAIJ_3_NaturalOrdering; 976 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_3_NaturalOrdering; 977 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=3\n");CHKERRQ(ierr); 978 break; 979 case 4: 980 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 981 inA->ops->solve = MatSolve_SeqSBAIJ_4_NaturalOrdering; 982 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_4_NaturalOrdering; 983 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=4\n");CHKERRQ(ierr); 984 break; 985 case 5: 986 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 987 inA->ops->solve = MatSolve_SeqSBAIJ_5_NaturalOrdering; 988 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_5_NaturalOrdering; 989 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=5\n");CHKERRQ(ierr); 990 break; 991 case 6: 992 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 993 inA->ops->solve = MatSolve_SeqSBAIJ_6_NaturalOrdering; 994 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_6_NaturalOrdering; 995 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=6\n");CHKERRQ(ierr); 996 break; 997 case 7: 998 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 999 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_7_NaturalOrdering; 1000 inA->ops->solve = MatSolve_SeqSBAIJ_7_NaturalOrdering; 1001 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=7\n");CHKERRQ(ierr); 1002 break; 1003 default: 1004 a->row = row; 1005 a->icol = col; 1006 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1007 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1008 1009 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1010 ierr = ISInvertPermutation(col,PETSC_DECIDE, &(a->icol));CHKERRQ(ierr); 1011 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1012 1013 if (!a->solve_work) { 1014 ierr = PetscMalloc((A->rmap.N+a->rmap.bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1015 ierr = PetscLogObjectMemory(inA,(A->rmap.N+a->rmap.bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1016 } 1017 } 1018 1019 ierr = MatCholeskyFactorNumeric(inA,info,&outA);CHKERRQ(ierr); 1020 PetscFunctionReturn(0); 1021 } 1022 #endif 1023 1024 EXTERN_C_BEGIN 1025 #undef __FUNCT__ 1026 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1027 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1028 { 1029 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1030 PetscInt i,nz,n; 1031 1032 PetscFunctionBegin; 1033 nz = baij->maxnz; 1034 n = mat->cmap.n; 1035 for (i=0; i<nz; i++) { 1036 baij->j[i] = indices[i]; 1037 } 1038 baij->nz = nz; 1039 for (i=0; i<n; i++) { 1040 baij->ilen[i] = baij->imax[i]; 1041 } 1042 PetscFunctionReturn(0); 1043 } 1044 EXTERN_C_END 1045 1046 #undef __FUNCT__ 1047 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1048 /*@ 1049 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1050 in the matrix. 1051 1052 Input Parameters: 1053 + mat - the SeqSBAIJ matrix 1054 - indices - the column indices 1055 1056 Level: advanced 1057 1058 Notes: 1059 This can be called if you have precomputed the nonzero structure of the 1060 matrix and want to provide it to the matrix object to improve the performance 1061 of the MatSetValues() operation. 1062 1063 You MUST have set the correct numbers of nonzeros per row in the call to 1064 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1065 1066 MUST be called before any calls to MatSetValues() 1067 1068 .seealso: MatCreateSeqSBAIJ 1069 @*/ 1070 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1071 { 1072 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 1073 1074 PetscFunctionBegin; 1075 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 1076 PetscValidPointer(indices,2); 1077 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 1078 if (f) { 1079 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1080 } else { 1081 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 1082 } 1083 PetscFunctionReturn(0); 1084 } 1085 1086 #undef __FUNCT__ 1087 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1088 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1089 { 1090 PetscErrorCode ierr; 1091 1092 PetscFunctionBegin; 1093 /* If the two matrices have the same copy implementation, use fast copy. */ 1094 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1095 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1096 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1097 1098 if (a->i[A->rmap.N] != b->i[B->rmap.N]) { 1099 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1100 } 1101 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap.N])*sizeof(PetscScalar));CHKERRQ(ierr); 1102 } else { 1103 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1104 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1105 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1106 } 1107 PetscFunctionReturn(0); 1108 } 1109 1110 #undef __FUNCT__ 1111 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1112 PetscErrorCode MatSetUpPreallocation_SeqSBAIJ(Mat A) 1113 { 1114 PetscErrorCode ierr; 1115 1116 PetscFunctionBegin; 1117 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,PetscMax(A->rmap.bs,1),PETSC_DEFAULT,0);CHKERRQ(ierr); 1118 PetscFunctionReturn(0); 1119 } 1120 1121 #undef __FUNCT__ 1122 #define __FUNCT__ "MatGetArray_SeqSBAIJ" 1123 PetscErrorCode MatGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1124 { 1125 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1126 PetscFunctionBegin; 1127 *array = a->a; 1128 PetscFunctionReturn(0); 1129 } 1130 1131 #undef __FUNCT__ 1132 #define __FUNCT__ "MatRestoreArray_SeqSBAIJ" 1133 PetscErrorCode MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1134 { 1135 PetscFunctionBegin; 1136 PetscFunctionReturn(0); 1137 } 1138 1139 #include "petscblaslapack.h" 1140 #undef __FUNCT__ 1141 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1142 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1143 { 1144 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ *)X->data, *y=(Mat_SeqSBAIJ *)Y->data; 1145 PetscErrorCode ierr; 1146 PetscInt i,bs=Y->rmap.bs,bs2,j; 1147 PetscBLASInt bnz = (PetscBLASInt)x->nz,one = 1; 1148 1149 PetscFunctionBegin; 1150 if (str == SAME_NONZERO_PATTERN) { 1151 PetscScalar alpha = a; 1152 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 1153 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1154 if (y->xtoy && y->XtoY != X) { 1155 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 1156 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 1157 } 1158 if (!y->xtoy) { /* get xtoy */ 1159 ierr = MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 1160 y->XtoY = X; 1161 } 1162 bs2 = bs*bs; 1163 for (i=0; i<x->nz; i++) { 1164 j = 0; 1165 while (j < bs2){ 1166 y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]); 1167 j++; 1168 } 1169 } 1170 ierr = PetscInfo3(0,"ratio of nnz_s(X)/nnz_s(Y): %D/%D = %G\n",bs2*x->nz,bs2*y->nz,(PetscReal)(bs2*x->nz)/(bs2*y->nz));CHKERRQ(ierr); 1171 } else { 1172 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1173 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1174 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1175 } 1176 PetscFunctionReturn(0); 1177 } 1178 1179 #undef __FUNCT__ 1180 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1181 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg) 1182 { 1183 PetscFunctionBegin; 1184 *flg = PETSC_TRUE; 1185 PetscFunctionReturn(0); 1186 } 1187 1188 #undef __FUNCT__ 1189 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1190 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscTruth *flg) 1191 { 1192 PetscFunctionBegin; 1193 *flg = PETSC_TRUE; 1194 PetscFunctionReturn(0); 1195 } 1196 1197 #undef __FUNCT__ 1198 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1199 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscTruth *flg) 1200 { 1201 PetscFunctionBegin; 1202 *flg = PETSC_FALSE; 1203 PetscFunctionReturn(0); 1204 } 1205 1206 #undef __FUNCT__ 1207 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1208 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1209 { 1210 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1211 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1212 PetscScalar *aa = a->a; 1213 1214 PetscFunctionBegin; 1215 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1216 PetscFunctionReturn(0); 1217 } 1218 1219 #undef __FUNCT__ 1220 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1221 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1222 { 1223 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1224 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1225 PetscScalar *aa = a->a; 1226 1227 PetscFunctionBegin; 1228 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1229 PetscFunctionReturn(0); 1230 } 1231 1232 /* -------------------------------------------------------------------*/ 1233 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1234 MatGetRow_SeqSBAIJ, 1235 MatRestoreRow_SeqSBAIJ, 1236 MatMult_SeqSBAIJ_N, 1237 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1238 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1239 MatMultAdd_SeqSBAIJ_N, 1240 MatSolve_SeqSBAIJ_N, 1241 0, 1242 0, 1243 /*10*/ 0, 1244 0, 1245 MatCholeskyFactor_SeqSBAIJ, 1246 MatRelax_SeqSBAIJ, 1247 MatTranspose_SeqSBAIJ, 1248 /*15*/ MatGetInfo_SeqSBAIJ, 1249 MatEqual_SeqSBAIJ, 1250 MatGetDiagonal_SeqSBAIJ, 1251 MatDiagonalScale_SeqSBAIJ, 1252 MatNorm_SeqSBAIJ, 1253 /*20*/ 0, 1254 MatAssemblyEnd_SeqSBAIJ, 1255 0, 1256 MatSetOption_SeqSBAIJ, 1257 MatZeroEntries_SeqSBAIJ, 1258 /*25*/ 0, 1259 0, 1260 0, 1261 MatCholeskyFactorSymbolic_SeqSBAIJ, 1262 MatCholeskyFactorNumeric_SeqSBAIJ_N, 1263 /*30*/ MatSetUpPreallocation_SeqSBAIJ, 1264 0, 1265 MatICCFactorSymbolic_SeqSBAIJ, 1266 MatGetArray_SeqSBAIJ, 1267 MatRestoreArray_SeqSBAIJ, 1268 /*35*/ MatDuplicate_SeqSBAIJ, 1269 0, 1270 0, 1271 0, 1272 0, 1273 /*40*/ MatAXPY_SeqSBAIJ, 1274 MatGetSubMatrices_SeqSBAIJ, 1275 MatIncreaseOverlap_SeqSBAIJ, 1276 MatGetValues_SeqSBAIJ, 1277 MatCopy_SeqSBAIJ, 1278 /*45*/ 0, 1279 MatScale_SeqSBAIJ, 1280 0, 1281 0, 1282 0, 1283 /*50*/ 0, 1284 MatGetRowIJ_SeqSBAIJ, 1285 MatRestoreRowIJ_SeqSBAIJ, 1286 0, 1287 0, 1288 /*55*/ 0, 1289 0, 1290 0, 1291 0, 1292 MatSetValuesBlocked_SeqSBAIJ, 1293 /*60*/ MatGetSubMatrix_SeqSBAIJ, 1294 0, 1295 0, 1296 0, 1297 0, 1298 /*65*/ 0, 1299 0, 1300 0, 1301 0, 1302 0, 1303 /*70*/ MatGetRowMax_SeqSBAIJ, 1304 0, 1305 0, 1306 0, 1307 0, 1308 /*75*/ 0, 1309 0, 1310 0, 1311 0, 1312 0, 1313 /*80*/ 0, 1314 0, 1315 0, 1316 #if !defined(PETSC_USE_COMPLEX) 1317 MatGetInertia_SeqSBAIJ, 1318 #else 1319 0, 1320 #endif 1321 MatLoad_SeqSBAIJ, 1322 /*85*/ MatIsSymmetric_SeqSBAIJ, 1323 MatIsHermitian_SeqSBAIJ, 1324 MatIsStructurallySymmetric_SeqSBAIJ, 1325 0, 1326 0, 1327 /*90*/ 0, 1328 0, 1329 0, 1330 0, 1331 0, 1332 /*95*/ 0, 1333 0, 1334 0, 1335 0, 1336 0, 1337 /*100*/0, 1338 0, 1339 0, 1340 0, 1341 0, 1342 /*105*/0, 1343 MatRealPart_SeqSBAIJ, 1344 MatImaginaryPart_SeqSBAIJ, 1345 MatGetRowUpperTriangular_SeqSBAIJ, 1346 MatRestoreRowUpperTriangular_SeqSBAIJ 1347 }; 1348 1349 EXTERN_C_BEGIN 1350 #undef __FUNCT__ 1351 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1352 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqSBAIJ(Mat mat) 1353 { 1354 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1355 PetscInt nz = aij->i[mat->rmap.N]*mat->rmap.bs*aij->bs2; 1356 PetscErrorCode ierr; 1357 1358 PetscFunctionBegin; 1359 if (aij->nonew != 1) { 1360 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1361 } 1362 1363 /* allocate space for values if not already there */ 1364 if (!aij->saved_values) { 1365 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1366 } 1367 1368 /* copy values over */ 1369 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1370 PetscFunctionReturn(0); 1371 } 1372 EXTERN_C_END 1373 1374 EXTERN_C_BEGIN 1375 #undef __FUNCT__ 1376 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1377 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqSBAIJ(Mat mat) 1378 { 1379 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1380 PetscErrorCode ierr; 1381 PetscInt nz = aij->i[mat->rmap.N]*mat->rmap.bs*aij->bs2; 1382 1383 PetscFunctionBegin; 1384 if (aij->nonew != 1) { 1385 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1386 } 1387 if (!aij->saved_values) { 1388 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1389 } 1390 1391 /* copy values over */ 1392 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1393 PetscFunctionReturn(0); 1394 } 1395 EXTERN_C_END 1396 1397 EXTERN_C_BEGIN 1398 #undef __FUNCT__ 1399 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1400 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1401 { 1402 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1403 PetscErrorCode ierr; 1404 PetscInt i,mbs,bs2; 1405 PetscTruth skipallocation = PETSC_FALSE,flg; 1406 1407 PetscFunctionBegin; 1408 B->preallocated = PETSC_TRUE; 1409 ierr = PetscOptionsGetInt(B->prefix,"-mat_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1410 B->rmap.bs = B->cmap.bs = bs; 1411 ierr = PetscMapInitialize(B->comm,&B->rmap);CHKERRQ(ierr); 1412 ierr = PetscMapInitialize(B->comm,&B->cmap);CHKERRQ(ierr); 1413 1414 mbs = B->rmap.N/bs; 1415 bs2 = bs*bs; 1416 1417 if (mbs*bs != B->rmap.N) { 1418 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1419 } 1420 1421 if (nz == MAT_SKIP_ALLOCATION) { 1422 skipallocation = PETSC_TRUE; 1423 nz = 0; 1424 } 1425 1426 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1427 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1428 if (nnz) { 1429 for (i=0; i<mbs; i++) { 1430 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]); 1431 if (nnz[i] > mbs) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %D value %D rowlength %D",i,nnz[i],mbs); 1432 } 1433 } 1434 1435 ierr = PetscOptionsHasName(B->prefix,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1436 if (!flg) { 1437 switch (bs) { 1438 case 1: 1439 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1440 B->ops->solve = MatSolve_SeqSBAIJ_1; 1441 B->ops->solves = MatSolves_SeqSBAIJ_1; 1442 B->ops->solvetranspose = MatSolve_SeqSBAIJ_1; 1443 B->ops->mult = MatMult_SeqSBAIJ_1; 1444 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1445 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1446 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1447 break; 1448 case 2: 1449 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1450 B->ops->solve = MatSolve_SeqSBAIJ_2; 1451 B->ops->solvetranspose = MatSolve_SeqSBAIJ_2; 1452 B->ops->mult = MatMult_SeqSBAIJ_2; 1453 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1454 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1455 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1456 break; 1457 case 3: 1458 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1459 B->ops->solve = MatSolve_SeqSBAIJ_3; 1460 B->ops->solvetranspose = MatSolve_SeqSBAIJ_3; 1461 B->ops->mult = MatMult_SeqSBAIJ_3; 1462 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1463 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1464 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1465 break; 1466 case 4: 1467 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1468 B->ops->solve = MatSolve_SeqSBAIJ_4; 1469 B->ops->solvetranspose = MatSolve_SeqSBAIJ_4; 1470 B->ops->mult = MatMult_SeqSBAIJ_4; 1471 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1472 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1473 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1474 break; 1475 case 5: 1476 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1477 B->ops->solve = MatSolve_SeqSBAIJ_5; 1478 B->ops->solvetranspose = MatSolve_SeqSBAIJ_5; 1479 B->ops->mult = MatMult_SeqSBAIJ_5; 1480 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1481 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1482 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1483 break; 1484 case 6: 1485 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1486 B->ops->solve = MatSolve_SeqSBAIJ_6; 1487 B->ops->solvetranspose = MatSolve_SeqSBAIJ_6; 1488 B->ops->mult = MatMult_SeqSBAIJ_6; 1489 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1490 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1491 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1492 break; 1493 case 7: 1494 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1495 B->ops->solve = MatSolve_SeqSBAIJ_7; 1496 B->ops->solvetranspose = MatSolve_SeqSBAIJ_7; 1497 B->ops->mult = MatMult_SeqSBAIJ_7; 1498 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1499 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1500 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1501 break; 1502 } 1503 } 1504 1505 b->mbs = mbs; 1506 b->nbs = mbs; 1507 if (!skipallocation) { 1508 /* b->ilen will count nonzeros in each block row so far. */ 1509 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1510 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1511 if (!nnz) { 1512 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1513 else if (nz <= 0) nz = 1; 1514 for (i=0; i<mbs; i++) { 1515 b->imax[i] = nz; 1516 } 1517 nz = nz*mbs; /* total nz */ 1518 } else { 1519 nz = 0; 1520 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1521 } 1522 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1523 1524 /* allocate the matrix space */ 1525 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap.N+1,PetscInt,&b->i);CHKERRQ(ierr); 1526 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1527 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1528 b->singlemalloc = PETSC_TRUE; 1529 1530 /* pointer to beginning of each row */ 1531 b->i[0] = 0; 1532 for (i=1; i<mbs+1; i++) { 1533 b->i[i] = b->i[i-1] + b->imax[i-1]; 1534 } 1535 b->free_a = PETSC_TRUE; 1536 b->free_ij = PETSC_TRUE; 1537 } else { 1538 b->free_a = PETSC_FALSE; 1539 b->free_ij = PETSC_FALSE; 1540 } 1541 1542 B->rmap.bs = bs; 1543 b->bs2 = bs2; 1544 b->nz = 0; 1545 b->maxnz = nz*bs2; 1546 1547 b->inew = 0; 1548 b->jnew = 0; 1549 b->anew = 0; 1550 b->a2anew = 0; 1551 b->permute = PETSC_FALSE; 1552 PetscFunctionReturn(0); 1553 } 1554 EXTERN_C_END 1555 1556 EXTERN_C_BEGIN 1557 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1558 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1559 EXTERN_C_END 1560 1561 /*MC 1562 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1563 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1564 1565 Options Database Keys: 1566 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1567 1568 Level: beginner 1569 1570 .seealso: MatCreateSeqSBAIJ 1571 M*/ 1572 1573 EXTERN_C_BEGIN 1574 #undef __FUNCT__ 1575 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1576 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqSBAIJ(Mat B) 1577 { 1578 Mat_SeqSBAIJ *b; 1579 PetscErrorCode ierr; 1580 PetscMPIInt size; 1581 PetscTruth flg; 1582 1583 PetscFunctionBegin; 1584 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1585 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1586 1587 ierr = PetscNew(Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1588 B->data = (void*)b; 1589 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1590 B->ops->destroy = MatDestroy_SeqSBAIJ; 1591 B->ops->view = MatView_SeqSBAIJ; 1592 B->factor = 0; 1593 B->mapping = 0; 1594 b->row = 0; 1595 b->icol = 0; 1596 b->reallocs = 0; 1597 b->saved_values = 0; 1598 1599 1600 b->sorted = PETSC_FALSE; 1601 b->roworiented = PETSC_TRUE; 1602 b->nonew = 0; 1603 b->diag = 0; 1604 b->solve_work = 0; 1605 b->mult_work = 0; 1606 B->spptr = 0; 1607 b->keepzeroedrows = PETSC_FALSE; 1608 b->xtoy = 0; 1609 b->XtoY = 0; 1610 1611 b->inew = 0; 1612 b->jnew = 0; 1613 b->anew = 0; 1614 b->a2anew = 0; 1615 b->permute = PETSC_FALSE; 1616 1617 b->ignore_ltriangular = PETSC_FALSE; 1618 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_ignore_lower_triangular",&flg);CHKERRQ(ierr); 1619 if (flg) b->ignore_ltriangular = PETSC_TRUE; 1620 1621 b->getrow_utriangular = PETSC_FALSE; 1622 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_getrow_uppertriangular",&flg);CHKERRQ(ierr); 1623 if (flg) b->getrow_utriangular = PETSC_TRUE; 1624 1625 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1626 "MatStoreValues_SeqSBAIJ", 1627 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1628 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1629 "MatRetrieveValues_SeqSBAIJ", 1630 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1631 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1632 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1633 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1634 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1635 "MatConvert_SeqSBAIJ_SeqAIJ", 1636 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1637 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1638 "MatConvert_SeqSBAIJ_SeqBAIJ", 1639 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1640 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1641 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1642 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1643 1644 B->symmetric = PETSC_TRUE; 1645 B->structurally_symmetric = PETSC_TRUE; 1646 B->symmetric_set = PETSC_TRUE; 1647 B->structurally_symmetric_set = PETSC_TRUE; 1648 PetscFunctionReturn(0); 1649 } 1650 EXTERN_C_END 1651 1652 #undef __FUNCT__ 1653 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1654 /*@C 1655 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1656 compressed row) format. For good matrix assembly performance the 1657 user should preallocate the matrix storage by setting the parameter nz 1658 (or the array nnz). By setting these parameters accurately, performance 1659 during matrix assembly can be increased by more than a factor of 50. 1660 1661 Collective on Mat 1662 1663 Input Parameters: 1664 + A - the symmetric matrix 1665 . bs - size of block 1666 . nz - number of block nonzeros per block row (same for all rows) 1667 - nnz - array containing the number of block nonzeros in the upper triangular plus 1668 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1669 1670 Options Database Keys: 1671 . -mat_no_unroll - uses code that does not unroll the loops in the 1672 block calculations (much slower) 1673 . -mat_block_size - size of the blocks to use 1674 1675 Level: intermediate 1676 1677 Notes: 1678 Specify the preallocated storage with either nz or nnz (not both). 1679 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1680 allocation. For additional details, see the users manual chapter on 1681 matrices. 1682 1683 If the nnz parameter is given then the nz parameter is ignored 1684 1685 1686 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1687 @*/ 1688 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1689 { 1690 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]); 1691 1692 PetscFunctionBegin; 1693 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1694 if (f) { 1695 ierr = (*f)(B,bs,nz,nnz);CHKERRQ(ierr); 1696 } 1697 PetscFunctionReturn(0); 1698 } 1699 1700 #undef __FUNCT__ 1701 #define __FUNCT__ "MatCreateSeqSBAIJ" 1702 /*@C 1703 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1704 compressed row) format. For good matrix assembly performance the 1705 user should preallocate the matrix storage by setting the parameter nz 1706 (or the array nnz). By setting these parameters accurately, performance 1707 during matrix assembly can be increased by more than a factor of 50. 1708 1709 Collective on MPI_Comm 1710 1711 Input Parameters: 1712 + comm - MPI communicator, set to PETSC_COMM_SELF 1713 . bs - size of block 1714 . m - number of rows, or number of columns 1715 . nz - number of block nonzeros per block row (same for all rows) 1716 - nnz - array containing the number of block nonzeros in the upper triangular plus 1717 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1718 1719 Output Parameter: 1720 . A - the symmetric matrix 1721 1722 Options Database Keys: 1723 . -mat_no_unroll - uses code that does not unroll the loops in the 1724 block calculations (much slower) 1725 . -mat_block_size - size of the blocks to use 1726 1727 Level: intermediate 1728 1729 Notes: 1730 1731 Specify the preallocated storage with either nz or nnz (not both). 1732 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1733 allocation. For additional details, see the users manual chapter on 1734 matrices. 1735 1736 If the nnz parameter is given then the nz parameter is ignored 1737 1738 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1739 @*/ 1740 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 1741 { 1742 PetscErrorCode ierr; 1743 1744 PetscFunctionBegin; 1745 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1746 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1747 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1748 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 1749 PetscFunctionReturn(0); 1750 } 1751 1752 #undef __FUNCT__ 1753 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1754 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1755 { 1756 Mat C; 1757 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1758 PetscErrorCode ierr; 1759 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1760 1761 PetscFunctionBegin; 1762 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1763 1764 *B = 0; 1765 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1766 ierr = MatSetSizes(C,A->rmap.N,A->cmap.n,A->rmap.N,A->cmap.n);CHKERRQ(ierr); 1767 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1768 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1769 c = (Mat_SeqSBAIJ*)C->data; 1770 1771 C->preallocated = PETSC_TRUE; 1772 C->factor = A->factor; 1773 c->row = 0; 1774 c->icol = 0; 1775 c->saved_values = 0; 1776 c->keepzeroedrows = a->keepzeroedrows; 1777 C->assembled = PETSC_TRUE; 1778 1779 ierr = PetscMapCopy(A->comm,&A->rmap,&C->rmap);CHKERRQ(ierr); 1780 ierr = PetscMapCopy(A->comm,&A->cmap,&C->cmap);CHKERRQ(ierr); 1781 c->bs2 = a->bs2; 1782 c->mbs = a->mbs; 1783 c->nbs = a->nbs; 1784 1785 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 1786 for (i=0; i<mbs; i++) { 1787 c->imax[i] = a->imax[i]; 1788 c->ilen[i] = a->ilen[i]; 1789 } 1790 1791 /* allocate the matrix space */ 1792 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 1793 c->singlemalloc = PETSC_TRUE; 1794 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1795 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 1796 if (mbs > 0) { 1797 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1798 if (cpvalues == MAT_COPY_VALUES) { 1799 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1800 } else { 1801 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1802 } 1803 } 1804 1805 c->sorted = a->sorted; 1806 c->roworiented = a->roworiented; 1807 c->nonew = a->nonew; 1808 1809 if (a->diag) { 1810 ierr = PetscMalloc((mbs+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 1811 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1812 for (i=0; i<mbs; i++) { 1813 c->diag[i] = a->diag[i]; 1814 } 1815 } else c->diag = 0; 1816 c->nz = a->nz; 1817 c->maxnz = a->maxnz; 1818 c->solve_work = 0; 1819 c->mult_work = 0; 1820 c->free_a = PETSC_TRUE; 1821 c->free_ij = PETSC_TRUE; 1822 *B = C; 1823 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1824 PetscFunctionReturn(0); 1825 } 1826 1827 #undef __FUNCT__ 1828 #define __FUNCT__ "MatLoad_SeqSBAIJ" 1829 PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, MatType type,Mat *A) 1830 { 1831 Mat_SeqSBAIJ *a; 1832 Mat B; 1833 PetscErrorCode ierr; 1834 int fd; 1835 PetscMPIInt size; 1836 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 1837 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 1838 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows; 1839 PetscInt *masked,nmask,tmp,bs2,ishift; 1840 PetscScalar *aa; 1841 MPI_Comm comm = ((PetscObject)viewer)->comm; 1842 1843 PetscFunctionBegin; 1844 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1845 bs2 = bs*bs; 1846 1847 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1848 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1849 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1850 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1851 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1852 M = header[1]; N = header[2]; nz = header[3]; 1853 1854 if (header[3] < 0) { 1855 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 1856 } 1857 1858 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1859 1860 /* 1861 This code adds extra rows to make sure the number of rows is 1862 divisible by the blocksize 1863 */ 1864 mbs = M/bs; 1865 extra_rows = bs - M + bs*(mbs); 1866 if (extra_rows == bs) extra_rows = 0; 1867 else mbs++; 1868 if (extra_rows) { 1869 ierr = PetscInfo(0,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 1870 } 1871 1872 /* read in row lengths */ 1873 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 1874 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1875 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1876 1877 /* read in column indices */ 1878 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 1879 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1880 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1881 1882 /* loop over row lengths determining block row lengths */ 1883 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 1884 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 1885 ierr = PetscMalloc(2*mbs*sizeof(PetscInt),&mask);CHKERRQ(ierr); 1886 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 1887 masked = mask + mbs; 1888 rowcount = 0; nzcount = 0; 1889 for (i=0; i<mbs; i++) { 1890 nmask = 0; 1891 for (j=0; j<bs; j++) { 1892 kmax = rowlengths[rowcount]; 1893 for (k=0; k<kmax; k++) { 1894 tmp = jj[nzcount++]/bs; /* block col. index */ 1895 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 1896 } 1897 rowcount++; 1898 } 1899 s_browlengths[i] += nmask; 1900 1901 /* zero out the mask elements we set */ 1902 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1903 } 1904 1905 /* create our matrix */ 1906 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 1907 ierr = MatSetSizes(B,M+extra_rows,N+extra_rows,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 1908 ierr = MatSetType(B,type);CHKERRQ(ierr); 1909 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,s_browlengths);CHKERRQ(ierr); 1910 a = (Mat_SeqSBAIJ*)B->data; 1911 1912 /* set matrix "i" values */ 1913 a->i[0] = 0; 1914 for (i=1; i<= mbs; i++) { 1915 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 1916 a->ilen[i-1] = s_browlengths[i-1]; 1917 } 1918 a->nz = a->i[mbs]; 1919 1920 /* read in nonzero values */ 1921 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 1922 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1923 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1924 1925 /* set "a" and "j" values into matrix */ 1926 nzcount = 0; jcount = 0; 1927 for (i=0; i<mbs; i++) { 1928 nzcountb = nzcount; 1929 nmask = 0; 1930 for (j=0; j<bs; j++) { 1931 kmax = rowlengths[i*bs+j]; 1932 for (k=0; k<kmax; k++) { 1933 tmp = jj[nzcount++]/bs; /* block col. index */ 1934 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 1935 } 1936 } 1937 /* sort the masked values */ 1938 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1939 1940 /* set "j" values into matrix */ 1941 maskcount = 1; 1942 for (j=0; j<nmask; j++) { 1943 a->j[jcount++] = masked[j]; 1944 mask[masked[j]] = maskcount++; 1945 } 1946 1947 /* set "a" values into matrix */ 1948 ishift = bs2*a->i[i]; 1949 for (j=0; j<bs; j++) { 1950 kmax = rowlengths[i*bs+j]; 1951 for (k=0; k<kmax; k++) { 1952 tmp = jj[nzcountb]/bs ; /* block col. index */ 1953 if (tmp >= i){ 1954 block = mask[tmp] - 1; 1955 point = jj[nzcountb] - bs*tmp; 1956 idx = ishift + bs2*block + j + bs*point; 1957 a->a[idx] = aa[nzcountb]; 1958 } 1959 nzcountb++; 1960 } 1961 } 1962 /* zero out the mask elements we set */ 1963 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1964 } 1965 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 1966 1967 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1968 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 1969 ierr = PetscFree(aa);CHKERRQ(ierr); 1970 ierr = PetscFree(jj);CHKERRQ(ierr); 1971 ierr = PetscFree(mask);CHKERRQ(ierr); 1972 1973 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1974 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1975 ierr = MatView_Private(B);CHKERRQ(ierr); 1976 *A = B; 1977 PetscFunctionReturn(0); 1978 } 1979 1980 #undef __FUNCT__ 1981 #define __FUNCT__ "MatRelax_SeqSBAIJ" 1982 PetscErrorCode MatRelax_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1983 { 1984 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1985 MatScalar *aa=a->a,*v,*v1; 1986 PetscScalar *x,*b,*t,sum,d; 1987 PetscErrorCode ierr; 1988 PetscInt m=a->mbs,bs=A->rmap.bs,*ai=a->i,*aj=a->j; 1989 PetscInt nz,nz1,*vj,*vj1,i; 1990 1991 PetscFunctionBegin; 1992 its = its*lits; 1993 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 1994 1995 if (bs > 1) 1996 SETERRQ(PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented"); 1997 1998 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1999 if (xx != bb) { 2000 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 2001 } else { 2002 b = x; 2003 } 2004 2005 if (!a->relax_work) { 2006 ierr = PetscMalloc(m*sizeof(PetscScalar),&a->relax_work);CHKERRQ(ierr); 2007 } 2008 t = a->relax_work; 2009 2010 if (flag & SOR_ZERO_INITIAL_GUESS) { 2011 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2012 for (i=0; i<m; i++) 2013 t[i] = b[i]; 2014 2015 for (i=0; i<m; i++){ 2016 d = *(aa + ai[i]); /* diag[i] */ 2017 v = aa + ai[i] + 1; 2018 vj = aj + ai[i] + 1; 2019 nz = ai[i+1] - ai[i] - 1; 2020 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2021 x[i] = omega*t[i]/d; 2022 while (nz--) t[*vj++] -= x[i]*(*v++); /* update rhs */ 2023 } 2024 } 2025 2026 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2027 if (!(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)){ 2028 t = b; 2029 } 2030 2031 for (i=m-1; i>=0; i--){ 2032 d = *(aa + ai[i]); 2033 v = aa + ai[i] + 1; 2034 vj = aj + ai[i] + 1; 2035 nz = ai[i+1] - ai[i] - 1; 2036 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2037 sum = t[i]; 2038 while (nz--) sum -= x[*vj++]*(*v++); 2039 x[i] = (1-omega)*x[i] + omega*sum/d; 2040 } 2041 t = a->relax_work; 2042 } 2043 its--; 2044 } 2045 2046 while (its--) { 2047 /* 2048 forward sweep: 2049 for i=0,...,m-1: 2050 sum[i] = (b[i] - U(i,:)x )/d[i]; 2051 x[i] = (1-omega)x[i] + omega*sum[i]; 2052 b = b - x[i]*U^T(i,:); 2053 2054 */ 2055 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2056 for (i=0; i<m; i++) 2057 t[i] = b[i]; 2058 2059 for (i=0; i<m; i++){ 2060 d = *(aa + ai[i]); /* diag[i] */ 2061 v = aa + ai[i] + 1; v1=v; 2062 vj = aj + ai[i] + 1; vj1=vj; 2063 nz = ai[i+1] - ai[i] - 1; nz1=nz; 2064 sum = t[i]; 2065 ierr = PetscLogFlops(4*nz-2);CHKERRQ(ierr); 2066 while (nz1--) sum -= (*v1++)*x[*vj1++]; 2067 x[i] = (1-omega)*x[i] + omega*sum/d; 2068 while (nz--) t[*vj++] -= x[i]*(*v++); 2069 } 2070 } 2071 2072 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2073 /* 2074 backward sweep: 2075 b = b - x[i]*U^T(i,:), i=0,...,n-2 2076 for i=m-1,...,0: 2077 sum[i] = (b[i] - U(i,:)x )/d[i]; 2078 x[i] = (1-omega)x[i] + omega*sum[i]; 2079 */ 2080 /* if there was a forward sweep done above then I thing the next two for loops are not needed */ 2081 for (i=0; i<m; i++) 2082 t[i] = b[i]; 2083 2084 for (i=0; i<m-1; i++){ /* update rhs */ 2085 v = aa + ai[i] + 1; 2086 vj = aj + ai[i] + 1; 2087 nz = ai[i+1] - ai[i] - 1; 2088 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2089 while (nz--) t[*vj++] -= x[i]*(*v++); 2090 } 2091 for (i=m-1; i>=0; i--){ 2092 d = *(aa + ai[i]); 2093 v = aa + ai[i] + 1; 2094 vj = aj + ai[i] + 1; 2095 nz = ai[i+1] - ai[i] - 1; 2096 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2097 sum = t[i]; 2098 while (nz--) sum -= x[*vj++]*(*v++); 2099 x[i] = (1-omega)*x[i] + omega*sum/d; 2100 } 2101 } 2102 } 2103 2104 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 2105 if (bb != xx) { 2106 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 2107 } 2108 PetscFunctionReturn(0); 2109 } 2110 2111 #undef __FUNCT__ 2112 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2113 /*@ 2114 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2115 (upper triangular entries in CSR format) provided by the user. 2116 2117 Collective on MPI_Comm 2118 2119 Input Parameters: 2120 + comm - must be an MPI communicator of size 1 2121 . bs - size of block 2122 . m - number of rows 2123 . n - number of columns 2124 . i - row indices 2125 . j - column indices 2126 - a - matrix values 2127 2128 Output Parameter: 2129 . mat - the matrix 2130 2131 Level: intermediate 2132 2133 Notes: 2134 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2135 once the matrix is destroyed 2136 2137 You cannot set new nonzero locations into this matrix, that will generate an error. 2138 2139 The i and j indices are 0 based 2140 2141 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2142 2143 @*/ 2144 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2145 { 2146 PetscErrorCode ierr; 2147 PetscInt ii; 2148 Mat_SeqSBAIJ *sbaij; 2149 2150 PetscFunctionBegin; 2151 if (bs != 1) SETERRQ1(PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2152 if (i[0]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2153 2154 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2155 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2156 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2157 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2158 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2159 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2160 2161 sbaij->i = i; 2162 sbaij->j = j; 2163 sbaij->a = a; 2164 sbaij->singlemalloc = PETSC_FALSE; 2165 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2166 sbaij->free_a = PETSC_FALSE; 2167 sbaij->free_ij = PETSC_FALSE; 2168 2169 for (ii=0; ii<m; ii++) { 2170 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2171 #if defined(PETSC_USE_DEBUG) 2172 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 2173 #endif 2174 } 2175 #if defined(PETSC_USE_DEBUG) 2176 for (ii=0; ii<sbaij->i[m]; ii++) { 2177 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2178 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2179 } 2180 #endif 2181 2182 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2183 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2184 PetscFunctionReturn(0); 2185 } 2186 2187 2188 2189 2190 2191