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) {ierr = ISDestroy(a->row);CHKERRQ(ierr);} 103 if (a->col){ierr = ISDestroy(a->col);CHKERRQ(ierr);} 104 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 105 ierr = PetscFree(a->diag);CHKERRQ(ierr); 106 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 107 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 108 ierr = PetscFree(a->relax_work);CHKERRQ(ierr); 109 ierr = PetscFree(a->solves_work);CHKERRQ(ierr); 110 ierr = PetscFree(a->mult_work);CHKERRQ(ierr); 111 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 112 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 113 114 ierr = PetscFree(a->inew);CHKERRQ(ierr); 115 ierr = PetscFree(a);CHKERRQ(ierr); 116 117 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 118 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 119 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 120 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C","",PETSC_NULL);CHKERRQ(ierr); 121 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 122 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 123 PetscFunctionReturn(0); 124 } 125 126 #undef __FUNCT__ 127 #define __FUNCT__ "MatSetOption_SeqSBAIJ" 128 PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op) 129 { 130 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 131 PetscErrorCode ierr; 132 133 PetscFunctionBegin; 134 switch (op) { 135 case MAT_ROW_ORIENTED: 136 a->roworiented = PETSC_TRUE; 137 break; 138 case MAT_COLUMN_ORIENTED: 139 a->roworiented = PETSC_FALSE; 140 break; 141 case MAT_COLUMNS_SORTED: 142 a->sorted = PETSC_TRUE; 143 break; 144 case MAT_COLUMNS_UNSORTED: 145 a->sorted = PETSC_FALSE; 146 break; 147 case MAT_KEEP_ZEROED_ROWS: 148 a->keepzeroedrows = PETSC_TRUE; 149 break; 150 case MAT_NO_NEW_NONZERO_LOCATIONS: 151 a->nonew = 1; 152 break; 153 case MAT_NEW_NONZERO_LOCATION_ERR: 154 a->nonew = -1; 155 break; 156 case MAT_NEW_NONZERO_ALLOCATION_ERR: 157 a->nonew = -2; 158 break; 159 case MAT_YES_NEW_NONZERO_LOCATIONS: 160 a->nonew = 0; 161 break; 162 case MAT_ROWS_SORTED: 163 case MAT_ROWS_UNSORTED: 164 case MAT_YES_NEW_DIAGONALS: 165 case MAT_IGNORE_OFF_PROC_ENTRIES: 166 case MAT_USE_HASH_TABLE: 167 ierr = PetscInfo1(A,"Option %d ignored\n",op);CHKERRQ(ierr); 168 break; 169 case MAT_NO_NEW_DIAGONALS: 170 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 171 case MAT_NOT_SYMMETRIC: 172 case MAT_NOT_STRUCTURALLY_SYMMETRIC: 173 case MAT_HERMITIAN: 174 SETERRQ(PETSC_ERR_SUP,"Matrix must be symmetric"); 175 case MAT_SYMMETRIC: 176 case MAT_STRUCTURALLY_SYMMETRIC: 177 case MAT_NOT_HERMITIAN: 178 case MAT_SYMMETRY_ETERNAL: 179 case MAT_NOT_SYMMETRY_ETERNAL: 180 case MAT_IGNORE_LOWER_TRIANGULAR: 181 a->ignore_ltriangular = PETSC_TRUE; 182 break; 183 case MAT_ERROR_LOWER_TRIANGULAR: 184 a->ignore_ltriangular = PETSC_FALSE; 185 break; 186 case MAT_GETROW_UPPERTRIANGULAR: 187 a->getrow_utriangular = PETSC_TRUE; 188 break; 189 default: 190 SETERRQ1(PETSC_ERR_SUP,"unknown option %d",op); 191 } 192 PetscFunctionReturn(0); 193 } 194 195 #undef __FUNCT__ 196 #define __FUNCT__ "MatGetRow_SeqSBAIJ" 197 PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **v) 198 { 199 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 200 PetscErrorCode ierr; 201 PetscInt itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*cols_i,bs2; 202 MatScalar *aa,*aa_i; 203 PetscScalar *v_i; 204 205 PetscFunctionBegin; 206 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()"); 207 /* Get the upper triangular part of the row */ 208 bs = A->rmap.bs; 209 ai = a->i; 210 aj = a->j; 211 aa = a->a; 212 bs2 = a->bs2; 213 214 if (row < 0 || row >= A->rmap.N) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Row %D out of range", row); 215 216 bn = row/bs; /* Block number */ 217 bp = row % bs; /* Block position */ 218 M = ai[bn+1] - ai[bn]; 219 *ncols = bs*M; 220 221 if (v) { 222 *v = 0; 223 if (*ncols) { 224 ierr = PetscMalloc((*ncols+row)*sizeof(PetscScalar),v);CHKERRQ(ierr); 225 for (i=0; i<M; i++) { /* for each block in the block row */ 226 v_i = *v + i*bs; 227 aa_i = aa + bs2*(ai[bn] + i); 228 for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];} 229 } 230 } 231 } 232 233 if (cols) { 234 *cols = 0; 235 if (*ncols) { 236 ierr = PetscMalloc((*ncols+row)*sizeof(PetscInt),cols);CHKERRQ(ierr); 237 for (i=0; i<M; i++) { /* for each block in the block row */ 238 cols_i = *cols + i*bs; 239 itmp = bs*aj[ai[bn] + i]; 240 for (j=0; j<bs; j++) {cols_i[j] = itmp++;} 241 } 242 } 243 } 244 245 /*search column A(0:row-1,row) (=A(row,0:row-1)). Could be expensive! */ 246 /* this segment is currently removed, so only entries in the upper triangle are obtained */ 247 #ifdef column_search 248 v_i = *v + M*bs; 249 cols_i = *cols + M*bs; 250 for (i=0; i<bn; i++){ /* for each block row */ 251 M = ai[i+1] - ai[i]; 252 for (j=0; j<M; j++){ 253 itmp = aj[ai[i] + j]; /* block column value */ 254 if (itmp == bn){ 255 aa_i = aa + bs2*(ai[i] + j) + bs*bp; 256 for (k=0; k<bs; k++) { 257 *cols_i++ = i*bs+k; 258 *v_i++ = aa_i[k]; 259 } 260 *ncols += bs; 261 break; 262 } 263 } 264 } 265 #endif 266 PetscFunctionReturn(0); 267 } 268 269 #undef __FUNCT__ 270 #define __FUNCT__ "MatRestoreRow_SeqSBAIJ" 271 PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 272 { 273 PetscErrorCode ierr; 274 275 PetscFunctionBegin; 276 if (idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);} 277 if (v) {ierr = PetscFree(*v);CHKERRQ(ierr);} 278 PetscFunctionReturn(0); 279 } 280 281 #undef __FUNCT__ 282 #define __FUNCT__ "MatGetRowUpperTriangular_SeqSBAIJ" 283 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A) 284 { 285 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 286 287 PetscFunctionBegin; 288 a->getrow_utriangular = PETSC_TRUE; 289 PetscFunctionReturn(0); 290 } 291 #undef __FUNCT__ 292 #define __FUNCT__ "MatRestoreRowUpperTriangular_SeqSBAIJ" 293 PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A) 294 { 295 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 296 297 PetscFunctionBegin; 298 a->getrow_utriangular = PETSC_FALSE; 299 PetscFunctionReturn(0); 300 } 301 302 #undef __FUNCT__ 303 #define __FUNCT__ "MatTranspose_SeqSBAIJ" 304 PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,Mat *B) 305 { 306 PetscErrorCode ierr; 307 PetscFunctionBegin; 308 ierr = MatDuplicate(A,MAT_COPY_VALUES,B);CHKERRQ(ierr); 309 PetscFunctionReturn(0); 310 } 311 312 #undef __FUNCT__ 313 #define __FUNCT__ "MatView_SeqSBAIJ_ASCII" 314 static PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 315 { 316 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 317 PetscErrorCode ierr; 318 PetscInt i,j,bs = A->rmap.bs,k,l,bs2=a->bs2; 319 const char *name; 320 PetscViewerFormat format; 321 322 PetscFunctionBegin; 323 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 324 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 325 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 326 ierr = PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);CHKERRQ(ierr); 327 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 328 SETERRQ(PETSC_ERR_SUP,"Matlab format not supported"); 329 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 330 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 331 for (i=0; i<a->mbs; i++) { 332 for (j=0; j<bs; j++) { 333 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 334 for (k=a->i[i]; k<a->i[i+1]; k++) { 335 for (l=0; l<bs; l++) { 336 #if defined(PETSC_USE_COMPLEX) 337 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 338 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 339 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 340 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 341 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 342 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 343 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 344 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 345 } 346 #else 347 if (a->a[bs2*k + l*bs + j] != 0.0) { 348 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 349 } 350 #endif 351 } 352 } 353 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 354 } 355 } 356 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 357 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 358 PetscFunctionReturn(0); 359 } else { 360 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 361 for (i=0; i<a->mbs; i++) { 362 for (j=0; j<bs; j++) { 363 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 364 for (k=a->i[i]; k<a->i[i+1]; k++) { 365 for (l=0; l<bs; l++) { 366 #if defined(PETSC_USE_COMPLEX) 367 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 368 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 369 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 370 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 371 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 372 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 373 } else { 374 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 375 } 376 #else 377 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 378 #endif 379 } 380 } 381 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 382 } 383 } 384 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 385 } 386 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 387 PetscFunctionReturn(0); 388 } 389 390 #undef __FUNCT__ 391 #define __FUNCT__ "MatView_SeqSBAIJ_Draw_Zoom" 392 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 393 { 394 Mat A = (Mat) Aa; 395 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 396 PetscErrorCode ierr; 397 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap.bs,bs2=a->bs2; 398 PetscMPIInt rank; 399 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 400 MatScalar *aa; 401 MPI_Comm comm; 402 PetscViewer viewer; 403 404 PetscFunctionBegin; 405 /* 406 This is nasty. If this is called from an originally parallel matrix 407 then all processes call this,but only the first has the matrix so the 408 rest should return immediately. 409 */ 410 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 411 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 412 if (rank) PetscFunctionReturn(0); 413 414 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 415 416 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 417 PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric"); 418 419 /* loop over matrix elements drawing boxes */ 420 color = PETSC_DRAW_BLUE; 421 for (i=0,row=0; i<mbs; i++,row+=bs) { 422 for (j=a->i[i]; j<a->i[i+1]; j++) { 423 y_l = A->rmap.N - row - 1.0; y_r = y_l + 1.0; 424 x_l = a->j[j]*bs; x_r = x_l + 1.0; 425 aa = a->a + j*bs2; 426 for (k=0; k<bs; k++) { 427 for (l=0; l<bs; l++) { 428 if (PetscRealPart(*aa++) >= 0.) continue; 429 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 430 } 431 } 432 } 433 } 434 color = PETSC_DRAW_CYAN; 435 for (i=0,row=0; i<mbs; i++,row+=bs) { 436 for (j=a->i[i]; j<a->i[i+1]; j++) { 437 y_l = A->rmap.N - row - 1.0; y_r = y_l + 1.0; 438 x_l = a->j[j]*bs; x_r = x_l + 1.0; 439 aa = a->a + j*bs2; 440 for (k=0; k<bs; k++) { 441 for (l=0; l<bs; l++) { 442 if (PetscRealPart(*aa++) != 0.) continue; 443 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 444 } 445 } 446 } 447 } 448 449 color = PETSC_DRAW_RED; 450 for (i=0,row=0; i<mbs; i++,row+=bs) { 451 for (j=a->i[i]; j<a->i[i+1]; j++) { 452 y_l = A->rmap.N - row - 1.0; y_r = y_l + 1.0; 453 x_l = a->j[j]*bs; x_r = x_l + 1.0; 454 aa = a->a + j*bs2; 455 for (k=0; k<bs; k++) { 456 for (l=0; l<bs; l++) { 457 if (PetscRealPart(*aa++) <= 0.) continue; 458 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 459 } 460 } 461 } 462 } 463 PetscFunctionReturn(0); 464 } 465 466 #undef __FUNCT__ 467 #define __FUNCT__ "MatView_SeqSBAIJ_Draw" 468 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 469 { 470 PetscErrorCode ierr; 471 PetscReal xl,yl,xr,yr,w,h; 472 PetscDraw draw; 473 PetscTruth isnull; 474 475 PetscFunctionBegin; 476 477 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 478 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 479 480 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 481 xr = A->rmap.N; yr = A->rmap.N; h = yr/10.0; w = xr/10.0; 482 xr += w; yr += h; xl = -w; yl = -h; 483 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 484 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 485 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 486 PetscFunctionReturn(0); 487 } 488 489 #undef __FUNCT__ 490 #define __FUNCT__ "MatView_SeqSBAIJ" 491 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 492 { 493 PetscErrorCode ierr; 494 PetscTruth iascii,isdraw; 495 496 PetscFunctionBegin; 497 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 498 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 499 if (iascii){ 500 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 501 } else if (isdraw) { 502 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 503 } else { 504 Mat B; 505 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 506 ierr = MatView(B,viewer);CHKERRQ(ierr); 507 ierr = MatDestroy(B);CHKERRQ(ierr); 508 } 509 PetscFunctionReturn(0); 510 } 511 512 513 #undef __FUNCT__ 514 #define __FUNCT__ "MatGetValues_SeqSBAIJ" 515 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 516 { 517 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 518 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 519 PetscInt *ai = a->i,*ailen = a->ilen; 520 PetscInt brow,bcol,ridx,cidx,bs=A->rmap.bs,bs2=a->bs2; 521 MatScalar *ap,*aa = a->a,zero = 0.0; 522 523 PetscFunctionBegin; 524 for (k=0; k<m; k++) { /* loop over rows */ 525 row = im[k]; brow = row/bs; 526 if (row < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); 527 if (row >= A->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap.N-1); 528 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 529 nrow = ailen[brow]; 530 for (l=0; l<n; l++) { /* loop over columns */ 531 if (in[l] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); 532 if (in[l] >= A->cmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap.n-1); 533 col = in[l] ; 534 bcol = col/bs; 535 cidx = col%bs; 536 ridx = row%bs; 537 high = nrow; 538 low = 0; /* assume unsorted */ 539 while (high-low > 5) { 540 t = (low+high)/2; 541 if (rp[t] > bcol) high = t; 542 else low = t; 543 } 544 for (i=low; i<high; i++) { 545 if (rp[i] > bcol) break; 546 if (rp[i] == bcol) { 547 *v++ = ap[bs2*i+bs*cidx+ridx]; 548 goto finished; 549 } 550 } 551 *v++ = zero; 552 finished:; 553 } 554 } 555 PetscFunctionReturn(0); 556 } 557 558 559 #undef __FUNCT__ 560 #define __FUNCT__ "MatSetValuesBlocked_SeqSBAIJ" 561 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 562 { 563 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 564 PetscErrorCode ierr; 565 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 566 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen; 567 PetscInt *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap.bs,stepval; 568 PetscTruth roworiented=a->roworiented; 569 const MatScalar *value = v; 570 MatScalar *ap,*aa = a->a,*bap; 571 572 PetscFunctionBegin; 573 if (roworiented) { 574 stepval = (n-1)*bs; 575 } else { 576 stepval = (m-1)*bs; 577 } 578 for (k=0; k<m; k++) { /* loop over added rows */ 579 row = im[k]; 580 if (row < 0) continue; 581 #if defined(PETSC_USE_DEBUG) 582 if (row >= a->mbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,a->mbs-1); 583 #endif 584 rp = aj + ai[row]; 585 ap = aa + bs2*ai[row]; 586 rmax = imax[row]; 587 nrow = ailen[row]; 588 low = 0; 589 high = nrow; 590 for (l=0; l<n; l++) { /* loop over added columns */ 591 if (in[l] < 0) continue; 592 col = in[l]; 593 #if defined(PETSC_USE_DEBUG) 594 if (col >= a->nbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",col,a->nbs-1); 595 #endif 596 if (col < row) continue; /* ignore lower triangular block */ 597 if (roworiented) { 598 value = v + k*(stepval+bs)*bs + l*bs; 599 } else { 600 value = v + l*(stepval+bs)*bs + k*bs; 601 } 602 if (col <= lastcol) low = 0; else high = nrow; 603 lastcol = col; 604 while (high-low > 7) { 605 t = (low+high)/2; 606 if (rp[t] > col) high = t; 607 else low = t; 608 } 609 for (i=low; i<high; i++) { 610 if (rp[i] > col) break; 611 if (rp[i] == col) { 612 bap = ap + bs2*i; 613 if (roworiented) { 614 if (is == ADD_VALUES) { 615 for (ii=0; ii<bs; ii++,value+=stepval) { 616 for (jj=ii; jj<bs2; jj+=bs) { 617 bap[jj] += *value++; 618 } 619 } 620 } else { 621 for (ii=0; ii<bs; ii++,value+=stepval) { 622 for (jj=ii; jj<bs2; jj+=bs) { 623 bap[jj] = *value++; 624 } 625 } 626 } 627 } else { 628 if (is == ADD_VALUES) { 629 for (ii=0; ii<bs; ii++,value+=stepval) { 630 for (jj=0; jj<bs; jj++) { 631 *bap++ += *value++; 632 } 633 } 634 } else { 635 for (ii=0; ii<bs; ii++,value+=stepval) { 636 for (jj=0; jj<bs; jj++) { 637 *bap++ = *value++; 638 } 639 } 640 } 641 } 642 goto noinsert2; 643 } 644 } 645 if (nonew == 1) goto noinsert2; 646 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 647 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew); 648 N = nrow++ - 1; high++; 649 /* shift up all the later entries in this row */ 650 for (ii=N; ii>=i; ii--) { 651 rp[ii+1] = rp[ii]; 652 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 653 } 654 if (N >= i) { 655 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 656 } 657 rp[i] = col; 658 bap = ap + bs2*i; 659 if (roworiented) { 660 for (ii=0; ii<bs; ii++,value+=stepval) { 661 for (jj=ii; jj<bs2; jj+=bs) { 662 bap[jj] = *value++; 663 } 664 } 665 } else { 666 for (ii=0; ii<bs; ii++,value+=stepval) { 667 for (jj=0; jj<bs; jj++) { 668 *bap++ = *value++; 669 } 670 } 671 } 672 noinsert2:; 673 low = i; 674 } 675 ailen[row] = nrow; 676 } 677 PetscFunctionReturn(0); 678 } 679 680 #undef __FUNCT__ 681 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 682 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 683 { 684 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 685 PetscErrorCode ierr; 686 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 687 PetscInt m = A->rmap.N,*ip,N,*ailen = a->ilen; 688 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 689 MatScalar *aa = a->a,*ap; 690 691 PetscFunctionBegin; 692 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 693 694 if (m) rmax = ailen[0]; 695 for (i=1; i<mbs; i++) { 696 /* move each row back by the amount of empty slots (fshift) before it*/ 697 fshift += imax[i-1] - ailen[i-1]; 698 rmax = PetscMax(rmax,ailen[i]); 699 if (fshift) { 700 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 701 N = ailen[i]; 702 for (j=0; j<N; j++) { 703 ip[j-fshift] = ip[j]; 704 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 705 } 706 } 707 ai[i] = ai[i-1] + ailen[i-1]; 708 } 709 if (mbs) { 710 fshift += imax[mbs-1] - ailen[mbs-1]; 711 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 712 } 713 /* reset ilen and imax for each row */ 714 for (i=0; i<mbs; i++) { 715 ailen[i] = imax[i] = ai[i+1] - ai[i]; 716 } 717 a->nz = ai[mbs]; 718 719 /* diagonals may have moved, reset it */ 720 if (a->diag) { 721 ierr = PetscMemcpy(a->diag,ai,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 722 } 723 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); 724 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 725 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 726 a->reallocs = 0; 727 A->info.nz_unneeded = (PetscReal)fshift*bs2; 728 PetscFunctionReturn(0); 729 } 730 731 /* 732 This function returns an array of flags which indicate the locations of contiguous 733 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 734 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 735 Assume: sizes should be long enough to hold all the values. 736 */ 737 #undef __FUNCT__ 738 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 739 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 740 { 741 PetscInt i,j,k,row; 742 PetscTruth flg; 743 744 PetscFunctionBegin; 745 for (i=0,j=0; i<n; j++) { 746 row = idx[i]; 747 if (row%bs!=0) { /* Not the begining of a block */ 748 sizes[j] = 1; 749 i++; 750 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 751 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 752 i++; 753 } else { /* Begining of the block, so check if the complete block exists */ 754 flg = PETSC_TRUE; 755 for (k=1; k<bs; k++) { 756 if (row+k != idx[i+k]) { /* break in the block */ 757 flg = PETSC_FALSE; 758 break; 759 } 760 } 761 if (flg) { /* No break in the bs */ 762 sizes[j] = bs; 763 i+= bs; 764 } else { 765 sizes[j] = 1; 766 i++; 767 } 768 } 769 } 770 *bs_max = j; 771 PetscFunctionReturn(0); 772 } 773 774 775 /* Only add/insert a(i,j) with i<=j (blocks). 776 Any a(i,j) with i>j input by user is ingored. 777 */ 778 779 #undef __FUNCT__ 780 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 781 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 782 { 783 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 784 PetscErrorCode ierr; 785 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 786 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 787 PetscInt *aj=a->j,nonew=a->nonew,bs=A->rmap.bs,brow,bcol; 788 PetscInt ridx,cidx,bs2=a->bs2; 789 MatScalar *ap,value,*aa=a->a,*bap; 790 791 PetscFunctionBegin; 792 for (k=0; k<m; k++) { /* loop over added rows */ 793 row = im[k]; /* row number */ 794 brow = row/bs; /* block row number */ 795 if (row < 0) continue; 796 #if defined(PETSC_USE_DEBUG) 797 if (row >= A->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap.N-1); 798 #endif 799 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 800 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 801 rmax = imax[brow]; /* maximum space allocated for this row */ 802 nrow = ailen[brow]; /* actual length of this row */ 803 low = 0; 804 805 for (l=0; l<n; l++) { /* loop over added columns */ 806 if (in[l] < 0) continue; 807 #if defined(PETSC_USE_DEBUG) 808 if (in[l] >= A->rmap.N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->rmap.N-1); 809 #endif 810 col = in[l]; 811 bcol = col/bs; /* block col number */ 812 813 if (brow > bcol) { 814 if (a->ignore_ltriangular){ 815 continue; /* ignore lower triangular values */ 816 } else { 817 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)"); 818 } 819 } 820 821 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 822 if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ 823 /* element value a(k,l) */ 824 if (roworiented) { 825 value = v[l + k*n]; 826 } else { 827 value = v[k + l*m]; 828 } 829 830 /* move pointer bap to a(k,l) quickly and add/insert value */ 831 if (col <= lastcol) low = 0; high = nrow; 832 lastcol = col; 833 while (high-low > 7) { 834 t = (low+high)/2; 835 if (rp[t] > bcol) high = t; 836 else low = t; 837 } 838 for (i=low; i<high; i++) { 839 if (rp[i] > bcol) break; 840 if (rp[i] == bcol) { 841 bap = ap + bs2*i + bs*cidx + ridx; 842 if (is == ADD_VALUES) *bap += value; 843 else *bap = value; 844 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 845 if (brow == bcol && ridx < cidx){ 846 bap = ap + bs2*i + bs*ridx + cidx; 847 if (is == ADD_VALUES) *bap += value; 848 else *bap = value; 849 } 850 goto noinsert1; 851 } 852 } 853 854 if (nonew == 1) goto noinsert1; 855 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 856 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew); 857 858 N = nrow++ - 1; high++; 859 /* shift up all the later entries in this row */ 860 for (ii=N; ii>=i; ii--) { 861 rp[ii+1] = rp[ii]; 862 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 863 } 864 if (N>=i) { 865 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 866 } 867 rp[i] = bcol; 868 ap[bs2*i + bs*cidx + ridx] = value; 869 noinsert1:; 870 low = i; 871 } 872 } /* end of loop over added columns */ 873 ailen[brow] = nrow; 874 } /* end of loop over added rows */ 875 PetscFunctionReturn(0); 876 } 877 878 EXTERN PetscErrorCode MatCholeskyFactorSymbolic_SeqSBAIJ(Mat,IS,MatFactorInfo*,Mat*); 879 EXTERN PetscErrorCode MatCholeskyFactor_SeqSBAIJ(Mat,IS,MatFactorInfo*); 880 EXTERN PetscErrorCode MatIncreaseOverlap_SeqSBAIJ(Mat,PetscInt,IS[],PetscInt); 881 EXTERN PetscErrorCode MatGetSubMatrix_SeqSBAIJ(Mat,IS,IS,PetscInt,MatReuse,Mat*); 882 EXTERN PetscErrorCode MatGetSubMatrices_SeqSBAIJ(Mat,PetscInt,const IS[],const IS[],MatReuse,Mat*[]); 883 EXTERN PetscErrorCode MatScale_SeqSBAIJ(Mat,PetscScalar); 884 EXTERN PetscErrorCode MatNorm_SeqSBAIJ(Mat,NormType,PetscReal *); 885 EXTERN PetscErrorCode MatEqual_SeqSBAIJ(Mat,Mat,PetscTruth*); 886 EXTERN PetscErrorCode MatGetDiagonal_SeqSBAIJ(Mat,Vec); 887 EXTERN PetscErrorCode MatDiagonalScale_SeqSBAIJ(Mat,Vec,Vec); 888 EXTERN PetscErrorCode MatGetInfo_SeqSBAIJ(Mat,MatInfoType,MatInfo *); 889 EXTERN PetscErrorCode MatZeroEntries_SeqSBAIJ(Mat); 890 EXTERN PetscErrorCode MatGetRowMax_SeqSBAIJ(Mat,Vec); 891 892 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_N(Mat,Vec,Vec); 893 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_1(Mat,Vec,Vec); 894 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_2(Mat,Vec,Vec); 895 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_3(Mat,Vec,Vec); 896 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_4(Mat,Vec,Vec); 897 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_5(Mat,Vec,Vec); 898 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_6(Mat,Vec,Vec); 899 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_7(Mat,Vec,Vec); 900 901 EXTERN PetscErrorCode MatSolves_SeqSBAIJ_1(Mat,Vecs,Vecs); 902 903 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_1_NaturalOrdering(Mat,Vec,Vec); 904 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_2_NaturalOrdering(Mat,Vec,Vec); 905 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_3_NaturalOrdering(Mat,Vec,Vec); 906 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_4_NaturalOrdering(Mat,Vec,Vec); 907 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_5_NaturalOrdering(Mat,Vec,Vec); 908 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_6_NaturalOrdering(Mat,Vec,Vec); 909 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_7_NaturalOrdering(Mat,Vec,Vec); 910 EXTERN PetscErrorCode MatSolve_SeqSBAIJ_N_NaturalOrdering(Mat,Vec,Vec); 911 912 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat,MatFactorInfo*,Mat*); 913 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_1(Mat,MatFactorInfo*,Mat*); 914 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat,MatFactorInfo*,Mat*); 915 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_3(Mat,MatFactorInfo*,Mat*); 916 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_4(Mat,MatFactorInfo*,Mat*); 917 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_5(Mat,MatFactorInfo*,Mat*); 918 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_6(Mat,MatFactorInfo*,Mat*); 919 EXTERN PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_7(Mat,MatFactorInfo*,Mat*); 920 EXTERN PetscErrorCode MatGetInertia_SeqSBAIJ(Mat,PetscInt*,PetscInt*,PetscInt*); 921 922 EXTERN PetscErrorCode MatMult_SeqSBAIJ_1(Mat,Vec,Vec); 923 EXTERN PetscErrorCode MatMult_SeqSBAIJ_2(Mat,Vec,Vec); 924 EXTERN PetscErrorCode MatMult_SeqSBAIJ_3(Mat,Vec,Vec); 925 EXTERN PetscErrorCode MatMult_SeqSBAIJ_4(Mat,Vec,Vec); 926 EXTERN PetscErrorCode MatMult_SeqSBAIJ_5(Mat,Vec,Vec); 927 EXTERN PetscErrorCode MatMult_SeqSBAIJ_6(Mat,Vec,Vec); 928 EXTERN PetscErrorCode MatMult_SeqSBAIJ_7(Mat,Vec,Vec); 929 EXTERN PetscErrorCode MatMult_SeqSBAIJ_N(Mat,Vec,Vec); 930 931 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_1(Mat,Vec,Vec,Vec); 932 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_2(Mat,Vec,Vec,Vec); 933 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_3(Mat,Vec,Vec,Vec); 934 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_4(Mat,Vec,Vec,Vec); 935 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_5(Mat,Vec,Vec,Vec); 936 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_6(Mat,Vec,Vec,Vec); 937 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_7(Mat,Vec,Vec,Vec); 938 EXTERN PetscErrorCode MatMultAdd_SeqSBAIJ_N(Mat,Vec,Vec,Vec); 939 940 #ifdef HAVE_MatICCFactor 941 /* modified from MatILUFactor_SeqSBAIJ, needs further work! */ 942 #undef __FUNCT__ 943 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 944 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,MatFactorInfo *info) 945 { 946 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 947 Mat outA; 948 PetscErrorCode ierr; 949 PetscTruth row_identity,col_identity; 950 951 PetscFunctionBegin; 952 outA = inA; 953 inA->factor = FACTOR_CHOLESKY; 954 955 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 956 /* 957 Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver 958 for ILU(0) factorization with natural ordering 959 */ 960 switch (a->rmap.bs) { 961 case 1: 962 inA->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering; 963 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering; 964 inA->ops->solves = MatSolves_SeqSBAIJ_1; 965 ierr = PetscInfo((inA,"Using special in-place natural ordering solvetrans BS=1\n");CHKERRQ(ierr); 966 case 2: 967 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 968 inA->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering; 969 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering; 970 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=2\n");CHKERRQ(ierr); 971 break; 972 case 3: 973 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 974 inA->ops->solve = MatSolve_SeqSBAIJ_3_NaturalOrdering; 975 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_3_NaturalOrdering; 976 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=3\n");CHKERRQ(ierr); 977 break; 978 case 4: 979 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 980 inA->ops->solve = MatSolve_SeqSBAIJ_4_NaturalOrdering; 981 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_4_NaturalOrdering; 982 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=4\n");CHKERRQ(ierr); 983 break; 984 case 5: 985 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 986 inA->ops->solve = MatSolve_SeqSBAIJ_5_NaturalOrdering; 987 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_5_NaturalOrdering; 988 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=5\n");CHKERRQ(ierr); 989 break; 990 case 6: 991 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 992 inA->ops->solve = MatSolve_SeqSBAIJ_6_NaturalOrdering; 993 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_6_NaturalOrdering; 994 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=6\n");CHKERRQ(ierr); 995 break; 996 case 7: 997 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 998 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_7_NaturalOrdering; 999 inA->ops->solve = MatSolve_SeqSBAIJ_7_NaturalOrdering; 1000 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=7\n");CHKERRQ(ierr); 1001 break; 1002 default: 1003 a->row = row; 1004 a->icol = col; 1005 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1006 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1007 1008 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1009 ierr = ISInvertPermutation(col,PETSC_DECIDE, &(a->icol));CHKERRQ(ierr); 1010 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1011 1012 if (!a->solve_work) { 1013 ierr = PetscMalloc((A->rmap.N+a->rmap.bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1014 ierr = PetscLogObjectMemory(inA,(A->rmap.N+a->rmap.bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1015 } 1016 } 1017 1018 ierr = MatCholeskyFactorNumeric(inA,info,&outA);CHKERRQ(ierr); 1019 PetscFunctionReturn(0); 1020 } 1021 #endif 1022 1023 EXTERN_C_BEGIN 1024 #undef __FUNCT__ 1025 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1026 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1027 { 1028 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1029 PetscInt i,nz,n; 1030 1031 PetscFunctionBegin; 1032 nz = baij->maxnz; 1033 n = mat->cmap.n; 1034 for (i=0; i<nz; i++) { 1035 baij->j[i] = indices[i]; 1036 } 1037 baij->nz = nz; 1038 for (i=0; i<n; i++) { 1039 baij->ilen[i] = baij->imax[i]; 1040 } 1041 PetscFunctionReturn(0); 1042 } 1043 EXTERN_C_END 1044 1045 #undef __FUNCT__ 1046 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1047 /*@ 1048 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1049 in the matrix. 1050 1051 Input Parameters: 1052 + mat - the SeqSBAIJ matrix 1053 - indices - the column indices 1054 1055 Level: advanced 1056 1057 Notes: 1058 This can be called if you have precomputed the nonzero structure of the 1059 matrix and want to provide it to the matrix object to improve the performance 1060 of the MatSetValues() operation. 1061 1062 You MUST have set the correct numbers of nonzeros per row in the call to 1063 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1064 1065 MUST be called before any calls to MatSetValues() 1066 1067 .seealso: MatCreateSeqSBAIJ 1068 @*/ 1069 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1070 { 1071 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 1072 1073 PetscFunctionBegin; 1074 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 1075 PetscValidPointer(indices,2); 1076 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 1077 if (f) { 1078 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1079 } else { 1080 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 1081 } 1082 PetscFunctionReturn(0); 1083 } 1084 1085 #undef __FUNCT__ 1086 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1087 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1088 { 1089 PetscErrorCode ierr; 1090 1091 PetscFunctionBegin; 1092 /* If the two matrices have the same copy implementation, use fast copy. */ 1093 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1094 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1095 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1096 1097 if (a->i[A->rmap.N] != b->i[B->rmap.N]) { 1098 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1099 } 1100 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap.N])*sizeof(PetscScalar));CHKERRQ(ierr); 1101 } else { 1102 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1103 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1104 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1105 } 1106 PetscFunctionReturn(0); 1107 } 1108 1109 #undef __FUNCT__ 1110 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1111 PetscErrorCode MatSetUpPreallocation_SeqSBAIJ(Mat A) 1112 { 1113 PetscErrorCode ierr; 1114 1115 PetscFunctionBegin; 1116 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,PetscMax(A->rmap.bs,1),PETSC_DEFAULT,0);CHKERRQ(ierr); 1117 PetscFunctionReturn(0); 1118 } 1119 1120 #undef __FUNCT__ 1121 #define __FUNCT__ "MatGetArray_SeqSBAIJ" 1122 PetscErrorCode MatGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1123 { 1124 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1125 PetscFunctionBegin; 1126 *array = a->a; 1127 PetscFunctionReturn(0); 1128 } 1129 1130 #undef __FUNCT__ 1131 #define __FUNCT__ "MatRestoreArray_SeqSBAIJ" 1132 PetscErrorCode MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1133 { 1134 PetscFunctionBegin; 1135 PetscFunctionReturn(0); 1136 } 1137 1138 #include "petscblaslapack.h" 1139 #undef __FUNCT__ 1140 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1141 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1142 { 1143 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ *)X->data, *y=(Mat_SeqSBAIJ *)Y->data; 1144 PetscErrorCode ierr; 1145 PetscInt i,bs=Y->rmap.bs,bs2,j; 1146 PetscBLASInt bnz = (PetscBLASInt)x->nz,one = 1; 1147 1148 PetscFunctionBegin; 1149 if (str == SAME_NONZERO_PATTERN) { 1150 PetscScalar alpha = a; 1151 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 1152 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1153 if (y->xtoy && y->XtoY != X) { 1154 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 1155 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 1156 } 1157 if (!y->xtoy) { /* get xtoy */ 1158 ierr = MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 1159 y->XtoY = X; 1160 } 1161 bs2 = bs*bs; 1162 for (i=0; i<x->nz; i++) { 1163 j = 0; 1164 while (j < bs2){ 1165 y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]); 1166 j++; 1167 } 1168 } 1169 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); 1170 } else { 1171 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1172 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1173 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1174 } 1175 PetscFunctionReturn(0); 1176 } 1177 1178 #undef __FUNCT__ 1179 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1180 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg) 1181 { 1182 PetscFunctionBegin; 1183 *flg = PETSC_TRUE; 1184 PetscFunctionReturn(0); 1185 } 1186 1187 #undef __FUNCT__ 1188 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1189 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscTruth *flg) 1190 { 1191 PetscFunctionBegin; 1192 *flg = PETSC_TRUE; 1193 PetscFunctionReturn(0); 1194 } 1195 1196 #undef __FUNCT__ 1197 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1198 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscTruth *flg) 1199 { 1200 PetscFunctionBegin; 1201 *flg = PETSC_FALSE; 1202 PetscFunctionReturn(0); 1203 } 1204 1205 #undef __FUNCT__ 1206 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1207 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1208 { 1209 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1210 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1211 PetscScalar *aa = a->a; 1212 1213 PetscFunctionBegin; 1214 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1215 PetscFunctionReturn(0); 1216 } 1217 1218 #undef __FUNCT__ 1219 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1220 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1221 { 1222 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1223 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1224 PetscScalar *aa = a->a; 1225 1226 PetscFunctionBegin; 1227 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1228 PetscFunctionReturn(0); 1229 } 1230 1231 /* -------------------------------------------------------------------*/ 1232 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1233 MatGetRow_SeqSBAIJ, 1234 MatRestoreRow_SeqSBAIJ, 1235 MatMult_SeqSBAIJ_N, 1236 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1237 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1238 MatMultAdd_SeqSBAIJ_N, 1239 MatSolve_SeqSBAIJ_N, 1240 0, 1241 0, 1242 /*10*/ 0, 1243 0, 1244 MatCholeskyFactor_SeqSBAIJ, 1245 MatRelax_SeqSBAIJ, 1246 MatTranspose_SeqSBAIJ, 1247 /*15*/ MatGetInfo_SeqSBAIJ, 1248 MatEqual_SeqSBAIJ, 1249 MatGetDiagonal_SeqSBAIJ, 1250 MatDiagonalScale_SeqSBAIJ, 1251 MatNorm_SeqSBAIJ, 1252 /*20*/ 0, 1253 MatAssemblyEnd_SeqSBAIJ, 1254 0, 1255 MatSetOption_SeqSBAIJ, 1256 MatZeroEntries_SeqSBAIJ, 1257 /*25*/ 0, 1258 0, 1259 0, 1260 MatCholeskyFactorSymbolic_SeqSBAIJ, 1261 MatCholeskyFactorNumeric_SeqSBAIJ_N, 1262 /*30*/ MatSetUpPreallocation_SeqSBAIJ, 1263 0, 1264 MatICCFactorSymbolic_SeqSBAIJ, 1265 MatGetArray_SeqSBAIJ, 1266 MatRestoreArray_SeqSBAIJ, 1267 /*35*/ MatDuplicate_SeqSBAIJ, 1268 0, 1269 0, 1270 0, 1271 0, 1272 /*40*/ MatAXPY_SeqSBAIJ, 1273 MatGetSubMatrices_SeqSBAIJ, 1274 MatIncreaseOverlap_SeqSBAIJ, 1275 MatGetValues_SeqSBAIJ, 1276 MatCopy_SeqSBAIJ, 1277 /*45*/ 0, 1278 MatScale_SeqSBAIJ, 1279 0, 1280 0, 1281 0, 1282 /*50*/ 0, 1283 MatGetRowIJ_SeqSBAIJ, 1284 MatRestoreRowIJ_SeqSBAIJ, 1285 0, 1286 0, 1287 /*55*/ 0, 1288 0, 1289 0, 1290 0, 1291 MatSetValuesBlocked_SeqSBAIJ, 1292 /*60*/ MatGetSubMatrix_SeqSBAIJ, 1293 0, 1294 0, 1295 0, 1296 0, 1297 /*65*/ 0, 1298 0, 1299 0, 1300 0, 1301 0, 1302 /*70*/ MatGetRowMax_SeqSBAIJ, 1303 0, 1304 0, 1305 0, 1306 0, 1307 /*75*/ 0, 1308 0, 1309 0, 1310 0, 1311 0, 1312 /*80*/ 0, 1313 0, 1314 0, 1315 #if !defined(PETSC_USE_COMPLEX) 1316 MatGetInertia_SeqSBAIJ, 1317 #else 1318 0, 1319 #endif 1320 MatLoad_SeqSBAIJ, 1321 /*85*/ MatIsSymmetric_SeqSBAIJ, 1322 MatIsHermitian_SeqSBAIJ, 1323 MatIsStructurallySymmetric_SeqSBAIJ, 1324 0, 1325 0, 1326 /*90*/ 0, 1327 0, 1328 0, 1329 0, 1330 0, 1331 /*95*/ 0, 1332 0, 1333 0, 1334 0, 1335 0, 1336 /*100*/0, 1337 0, 1338 0, 1339 0, 1340 0, 1341 /*105*/0, 1342 MatRealPart_SeqSBAIJ, 1343 MatImaginaryPart_SeqSBAIJ, 1344 MatGetRowUpperTriangular_SeqSBAIJ, 1345 MatRestoreRowUpperTriangular_SeqSBAIJ 1346 }; 1347 1348 EXTERN_C_BEGIN 1349 #undef __FUNCT__ 1350 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1351 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqSBAIJ(Mat mat) 1352 { 1353 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1354 PetscInt nz = aij->i[mat->rmap.N]*mat->rmap.bs*aij->bs2; 1355 PetscErrorCode ierr; 1356 1357 PetscFunctionBegin; 1358 if (aij->nonew != 1) { 1359 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1360 } 1361 1362 /* allocate space for values if not already there */ 1363 if (!aij->saved_values) { 1364 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1365 } 1366 1367 /* copy values over */ 1368 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1369 PetscFunctionReturn(0); 1370 } 1371 EXTERN_C_END 1372 1373 EXTERN_C_BEGIN 1374 #undef __FUNCT__ 1375 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1376 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqSBAIJ(Mat mat) 1377 { 1378 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1379 PetscErrorCode ierr; 1380 PetscInt nz = aij->i[mat->rmap.N]*mat->rmap.bs*aij->bs2; 1381 1382 PetscFunctionBegin; 1383 if (aij->nonew != 1) { 1384 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1385 } 1386 if (!aij->saved_values) { 1387 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1388 } 1389 1390 /* copy values over */ 1391 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1392 PetscFunctionReturn(0); 1393 } 1394 EXTERN_C_END 1395 1396 EXTERN_C_BEGIN 1397 #undef __FUNCT__ 1398 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1399 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1400 { 1401 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1402 PetscErrorCode ierr; 1403 PetscInt i,mbs,bs2; 1404 PetscTruth skipallocation = PETSC_FALSE,flg; 1405 1406 PetscFunctionBegin; 1407 B->preallocated = PETSC_TRUE; 1408 ierr = PetscOptionsGetInt(B->prefix,"-mat_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1409 B->rmap.bs = B->cmap.bs = bs; 1410 ierr = PetscMapInitialize(B->comm,&B->rmap);CHKERRQ(ierr); 1411 ierr = PetscMapInitialize(B->comm,&B->cmap);CHKERRQ(ierr); 1412 1413 mbs = B->rmap.N/bs; 1414 bs2 = bs*bs; 1415 1416 if (mbs*bs != B->rmap.N) { 1417 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1418 } 1419 1420 if (nz == MAT_SKIP_ALLOCATION) { 1421 skipallocation = PETSC_TRUE; 1422 nz = 0; 1423 } 1424 1425 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1426 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1427 if (nnz) { 1428 for (i=0; i<mbs; i++) { 1429 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]); 1430 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); 1431 } 1432 } 1433 1434 ierr = PetscOptionsHasName(B->prefix,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1435 if (!flg) { 1436 switch (bs) { 1437 case 1: 1438 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1439 B->ops->solve = MatSolve_SeqSBAIJ_1; 1440 B->ops->solves = MatSolves_SeqSBAIJ_1; 1441 B->ops->solvetranspose = MatSolve_SeqSBAIJ_1; 1442 B->ops->mult = MatMult_SeqSBAIJ_1; 1443 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1444 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1445 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1446 break; 1447 case 2: 1448 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1449 B->ops->solve = MatSolve_SeqSBAIJ_2; 1450 B->ops->solvetranspose = MatSolve_SeqSBAIJ_2; 1451 B->ops->mult = MatMult_SeqSBAIJ_2; 1452 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1453 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1454 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1455 break; 1456 case 3: 1457 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1458 B->ops->solve = MatSolve_SeqSBAIJ_3; 1459 B->ops->solvetranspose = MatSolve_SeqSBAIJ_3; 1460 B->ops->mult = MatMult_SeqSBAIJ_3; 1461 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1462 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1463 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1464 break; 1465 case 4: 1466 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1467 B->ops->solve = MatSolve_SeqSBAIJ_4; 1468 B->ops->solvetranspose = MatSolve_SeqSBAIJ_4; 1469 B->ops->mult = MatMult_SeqSBAIJ_4; 1470 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1471 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1472 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1473 break; 1474 case 5: 1475 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1476 B->ops->solve = MatSolve_SeqSBAIJ_5; 1477 B->ops->solvetranspose = MatSolve_SeqSBAIJ_5; 1478 B->ops->mult = MatMult_SeqSBAIJ_5; 1479 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1480 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1481 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1482 break; 1483 case 6: 1484 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1485 B->ops->solve = MatSolve_SeqSBAIJ_6; 1486 B->ops->solvetranspose = MatSolve_SeqSBAIJ_6; 1487 B->ops->mult = MatMult_SeqSBAIJ_6; 1488 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1489 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1490 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1491 break; 1492 case 7: 1493 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1494 B->ops->solve = MatSolve_SeqSBAIJ_7; 1495 B->ops->solvetranspose = MatSolve_SeqSBAIJ_7; 1496 B->ops->mult = MatMult_SeqSBAIJ_7; 1497 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1498 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1499 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1500 break; 1501 } 1502 } 1503 1504 b->mbs = mbs; 1505 b->nbs = mbs; 1506 if (!skipallocation) { 1507 /* b->ilen will count nonzeros in each block row so far. */ 1508 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1509 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1510 if (!nnz) { 1511 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1512 else if (nz <= 0) nz = 1; 1513 for (i=0; i<mbs; i++) { 1514 b->imax[i] = nz; 1515 } 1516 nz = nz*mbs; /* total nz */ 1517 } else { 1518 nz = 0; 1519 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1520 } 1521 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1522 1523 /* allocate the matrix space */ 1524 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap.N+1,PetscInt,&b->i);CHKERRQ(ierr); 1525 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1526 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1527 b->singlemalloc = PETSC_TRUE; 1528 1529 /* pointer to beginning of each row */ 1530 b->i[0] = 0; 1531 for (i=1; i<mbs+1; i++) { 1532 b->i[i] = b->i[i-1] + b->imax[i-1]; 1533 } 1534 b->free_a = PETSC_TRUE; 1535 b->free_ij = PETSC_TRUE; 1536 } else { 1537 b->free_a = PETSC_FALSE; 1538 b->free_ij = PETSC_FALSE; 1539 } 1540 1541 B->rmap.bs = bs; 1542 b->bs2 = bs2; 1543 b->nz = 0; 1544 b->maxnz = nz*bs2; 1545 1546 b->inew = 0; 1547 b->jnew = 0; 1548 b->anew = 0; 1549 b->a2anew = 0; 1550 b->permute = PETSC_FALSE; 1551 PetscFunctionReturn(0); 1552 } 1553 EXTERN_C_END 1554 1555 EXTERN_C_BEGIN 1556 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1557 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1558 EXTERN_C_END 1559 1560 /*MC 1561 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1562 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1563 1564 Options Database Keys: 1565 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1566 1567 Level: beginner 1568 1569 .seealso: MatCreateSeqSBAIJ 1570 M*/ 1571 1572 EXTERN_C_BEGIN 1573 #undef __FUNCT__ 1574 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1575 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqSBAIJ(Mat B) 1576 { 1577 Mat_SeqSBAIJ *b; 1578 PetscErrorCode ierr; 1579 PetscMPIInt size; 1580 PetscTruth flg; 1581 1582 PetscFunctionBegin; 1583 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1584 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1585 1586 ierr = PetscNew(Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1587 B->data = (void*)b; 1588 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1589 B->ops->destroy = MatDestroy_SeqSBAIJ; 1590 B->ops->view = MatView_SeqSBAIJ; 1591 B->factor = 0; 1592 B->mapping = 0; 1593 b->row = 0; 1594 b->icol = 0; 1595 b->reallocs = 0; 1596 b->saved_values = 0; 1597 1598 1599 b->sorted = PETSC_FALSE; 1600 b->roworiented = PETSC_TRUE; 1601 b->nonew = 0; 1602 b->diag = 0; 1603 b->solve_work = 0; 1604 b->mult_work = 0; 1605 B->spptr = 0; 1606 b->keepzeroedrows = PETSC_FALSE; 1607 b->xtoy = 0; 1608 b->XtoY = 0; 1609 1610 b->inew = 0; 1611 b->jnew = 0; 1612 b->anew = 0; 1613 b->a2anew = 0; 1614 b->permute = PETSC_FALSE; 1615 1616 b->ignore_ltriangular = PETSC_FALSE; 1617 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_ignore_lower_triangular",&flg);CHKERRQ(ierr); 1618 if (flg) b->ignore_ltriangular = PETSC_TRUE; 1619 1620 b->getrow_utriangular = PETSC_FALSE; 1621 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_getrow_uppertriangular",&flg);CHKERRQ(ierr); 1622 if (flg) b->getrow_utriangular = PETSC_TRUE; 1623 1624 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1625 "MatStoreValues_SeqSBAIJ", 1626 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1627 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1628 "MatRetrieveValues_SeqSBAIJ", 1629 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1630 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1631 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1632 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1633 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1634 "MatConvert_SeqSBAIJ_SeqAIJ", 1635 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1636 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1637 "MatConvert_SeqSBAIJ_SeqBAIJ", 1638 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1639 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1640 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1641 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1642 1643 B->symmetric = PETSC_TRUE; 1644 B->structurally_symmetric = PETSC_TRUE; 1645 B->symmetric_set = PETSC_TRUE; 1646 B->structurally_symmetric_set = PETSC_TRUE; 1647 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 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