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 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1649 PetscFunctionReturn(0); 1650 } 1651 EXTERN_C_END 1652 1653 #undef __FUNCT__ 1654 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1655 /*@C 1656 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1657 compressed row) format. For good matrix assembly performance the 1658 user should preallocate the matrix storage by setting the parameter nz 1659 (or the array nnz). By setting these parameters accurately, performance 1660 during matrix assembly can be increased by more than a factor of 50. 1661 1662 Collective on Mat 1663 1664 Input Parameters: 1665 + A - the symmetric matrix 1666 . bs - size of block 1667 . nz - number of block nonzeros per block row (same for all rows) 1668 - nnz - array containing the number of block nonzeros in the upper triangular plus 1669 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1670 1671 Options Database Keys: 1672 . -mat_no_unroll - uses code that does not unroll the loops in the 1673 block calculations (much slower) 1674 . -mat_block_size - size of the blocks to use 1675 1676 Level: intermediate 1677 1678 Notes: 1679 Specify the preallocated storage with either nz or nnz (not both). 1680 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1681 allocation. For additional details, see the users manual chapter on 1682 matrices. 1683 1684 If the nnz parameter is given then the nz parameter is ignored 1685 1686 1687 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1688 @*/ 1689 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1690 { 1691 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]); 1692 1693 PetscFunctionBegin; 1694 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1695 if (f) { 1696 ierr = (*f)(B,bs,nz,nnz);CHKERRQ(ierr); 1697 } 1698 PetscFunctionReturn(0); 1699 } 1700 1701 #undef __FUNCT__ 1702 #define __FUNCT__ "MatCreateSeqSBAIJ" 1703 /*@C 1704 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1705 compressed row) format. For good matrix assembly performance the 1706 user should preallocate the matrix storage by setting the parameter nz 1707 (or the array nnz). By setting these parameters accurately, performance 1708 during matrix assembly can be increased by more than a factor of 50. 1709 1710 Collective on MPI_Comm 1711 1712 Input Parameters: 1713 + comm - MPI communicator, set to PETSC_COMM_SELF 1714 . bs - size of block 1715 . m - number of rows, or number of columns 1716 . nz - number of block nonzeros per block row (same for all rows) 1717 - nnz - array containing the number of block nonzeros in the upper triangular plus 1718 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1719 1720 Output Parameter: 1721 . A - the symmetric matrix 1722 1723 Options Database Keys: 1724 . -mat_no_unroll - uses code that does not unroll the loops in the 1725 block calculations (much slower) 1726 . -mat_block_size - size of the blocks to use 1727 1728 Level: intermediate 1729 1730 Notes: 1731 1732 Specify the preallocated storage with either nz or nnz (not both). 1733 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1734 allocation. For additional details, see the users manual chapter on 1735 matrices. 1736 1737 If the nnz parameter is given then the nz parameter is ignored 1738 1739 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1740 @*/ 1741 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 1742 { 1743 PetscErrorCode ierr; 1744 1745 PetscFunctionBegin; 1746 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1747 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1748 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1749 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 1750 PetscFunctionReturn(0); 1751 } 1752 1753 #undef __FUNCT__ 1754 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1755 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1756 { 1757 Mat C; 1758 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1759 PetscErrorCode ierr; 1760 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1761 1762 PetscFunctionBegin; 1763 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1764 1765 *B = 0; 1766 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1767 ierr = MatSetSizes(C,A->rmap.N,A->cmap.n,A->rmap.N,A->cmap.n);CHKERRQ(ierr); 1768 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1769 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1770 c = (Mat_SeqSBAIJ*)C->data; 1771 1772 C->preallocated = PETSC_TRUE; 1773 C->factor = A->factor; 1774 c->row = 0; 1775 c->icol = 0; 1776 c->saved_values = 0; 1777 c->keepzeroedrows = a->keepzeroedrows; 1778 C->assembled = PETSC_TRUE; 1779 1780 ierr = PetscMapCopy(A->comm,&A->rmap,&C->rmap);CHKERRQ(ierr); 1781 ierr = PetscMapCopy(A->comm,&A->cmap,&C->cmap);CHKERRQ(ierr); 1782 c->bs2 = a->bs2; 1783 c->mbs = a->mbs; 1784 c->nbs = a->nbs; 1785 1786 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 1787 for (i=0; i<mbs; i++) { 1788 c->imax[i] = a->imax[i]; 1789 c->ilen[i] = a->ilen[i]; 1790 } 1791 1792 /* allocate the matrix space */ 1793 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 1794 c->singlemalloc = PETSC_TRUE; 1795 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1796 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 1797 if (mbs > 0) { 1798 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1799 if (cpvalues == MAT_COPY_VALUES) { 1800 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1801 } else { 1802 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1803 } 1804 } 1805 1806 c->sorted = a->sorted; 1807 c->roworiented = a->roworiented; 1808 c->nonew = a->nonew; 1809 1810 if (a->diag) { 1811 ierr = PetscMalloc((mbs+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 1812 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1813 for (i=0; i<mbs; i++) { 1814 c->diag[i] = a->diag[i]; 1815 } 1816 } else c->diag = 0; 1817 c->nz = a->nz; 1818 c->maxnz = a->maxnz; 1819 c->solve_work = 0; 1820 c->mult_work = 0; 1821 c->free_a = PETSC_TRUE; 1822 c->free_ij = PETSC_TRUE; 1823 *B = C; 1824 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1825 PetscFunctionReturn(0); 1826 } 1827 1828 #undef __FUNCT__ 1829 #define __FUNCT__ "MatLoad_SeqSBAIJ" 1830 PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, MatType type,Mat *A) 1831 { 1832 Mat_SeqSBAIJ *a; 1833 Mat B; 1834 PetscErrorCode ierr; 1835 int fd; 1836 PetscMPIInt size; 1837 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 1838 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 1839 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows; 1840 PetscInt *masked,nmask,tmp,bs2,ishift; 1841 PetscScalar *aa; 1842 MPI_Comm comm = ((PetscObject)viewer)->comm; 1843 1844 PetscFunctionBegin; 1845 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1846 bs2 = bs*bs; 1847 1848 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1849 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1850 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1851 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1852 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1853 M = header[1]; N = header[2]; nz = header[3]; 1854 1855 if (header[3] < 0) { 1856 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 1857 } 1858 1859 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1860 1861 /* 1862 This code adds extra rows to make sure the number of rows is 1863 divisible by the blocksize 1864 */ 1865 mbs = M/bs; 1866 extra_rows = bs - M + bs*(mbs); 1867 if (extra_rows == bs) extra_rows = 0; 1868 else mbs++; 1869 if (extra_rows) { 1870 ierr = PetscInfo(0,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 1871 } 1872 1873 /* read in row lengths */ 1874 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 1875 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1876 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1877 1878 /* read in column indices */ 1879 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 1880 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1881 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1882 1883 /* loop over row lengths determining block row lengths */ 1884 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 1885 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 1886 ierr = PetscMalloc(2*mbs*sizeof(PetscInt),&mask);CHKERRQ(ierr); 1887 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 1888 masked = mask + mbs; 1889 rowcount = 0; nzcount = 0; 1890 for (i=0; i<mbs; i++) { 1891 nmask = 0; 1892 for (j=0; j<bs; j++) { 1893 kmax = rowlengths[rowcount]; 1894 for (k=0; k<kmax; k++) { 1895 tmp = jj[nzcount++]/bs; /* block col. index */ 1896 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 1897 } 1898 rowcount++; 1899 } 1900 s_browlengths[i] += nmask; 1901 1902 /* zero out the mask elements we set */ 1903 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1904 } 1905 1906 /* create our matrix */ 1907 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 1908 ierr = MatSetSizes(B,M+extra_rows,N+extra_rows,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 1909 ierr = MatSetType(B,type);CHKERRQ(ierr); 1910 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,s_browlengths);CHKERRQ(ierr); 1911 a = (Mat_SeqSBAIJ*)B->data; 1912 1913 /* set matrix "i" values */ 1914 a->i[0] = 0; 1915 for (i=1; i<= mbs; i++) { 1916 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 1917 a->ilen[i-1] = s_browlengths[i-1]; 1918 } 1919 a->nz = a->i[mbs]; 1920 1921 /* read in nonzero values */ 1922 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 1923 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1924 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1925 1926 /* set "a" and "j" values into matrix */ 1927 nzcount = 0; jcount = 0; 1928 for (i=0; i<mbs; i++) { 1929 nzcountb = nzcount; 1930 nmask = 0; 1931 for (j=0; j<bs; j++) { 1932 kmax = rowlengths[i*bs+j]; 1933 for (k=0; k<kmax; k++) { 1934 tmp = jj[nzcount++]/bs; /* block col. index */ 1935 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 1936 } 1937 } 1938 /* sort the masked values */ 1939 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1940 1941 /* set "j" values into matrix */ 1942 maskcount = 1; 1943 for (j=0; j<nmask; j++) { 1944 a->j[jcount++] = masked[j]; 1945 mask[masked[j]] = maskcount++; 1946 } 1947 1948 /* set "a" values into matrix */ 1949 ishift = bs2*a->i[i]; 1950 for (j=0; j<bs; j++) { 1951 kmax = rowlengths[i*bs+j]; 1952 for (k=0; k<kmax; k++) { 1953 tmp = jj[nzcountb]/bs ; /* block col. index */ 1954 if (tmp >= i){ 1955 block = mask[tmp] - 1; 1956 point = jj[nzcountb] - bs*tmp; 1957 idx = ishift + bs2*block + j + bs*point; 1958 a->a[idx] = aa[nzcountb]; 1959 } 1960 nzcountb++; 1961 } 1962 } 1963 /* zero out the mask elements we set */ 1964 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1965 } 1966 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 1967 1968 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1969 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 1970 ierr = PetscFree(aa);CHKERRQ(ierr); 1971 ierr = PetscFree(jj);CHKERRQ(ierr); 1972 ierr = PetscFree(mask);CHKERRQ(ierr); 1973 1974 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1975 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1976 ierr = MatView_Private(B);CHKERRQ(ierr); 1977 *A = B; 1978 PetscFunctionReturn(0); 1979 } 1980 1981 #undef __FUNCT__ 1982 #define __FUNCT__ "MatRelax_SeqSBAIJ" 1983 PetscErrorCode MatRelax_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1984 { 1985 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1986 MatScalar *aa=a->a,*v,*v1; 1987 PetscScalar *x,*b,*t,sum,d; 1988 PetscErrorCode ierr; 1989 PetscInt m=a->mbs,bs=A->rmap.bs,*ai=a->i,*aj=a->j; 1990 PetscInt nz,nz1,*vj,*vj1,i; 1991 1992 PetscFunctionBegin; 1993 its = its*lits; 1994 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 1995 1996 if (bs > 1) 1997 SETERRQ(PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented"); 1998 1999 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 2000 if (xx != bb) { 2001 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 2002 } else { 2003 b = x; 2004 } 2005 2006 if (!a->relax_work) { 2007 ierr = PetscMalloc(m*sizeof(PetscScalar),&a->relax_work);CHKERRQ(ierr); 2008 } 2009 t = a->relax_work; 2010 2011 if (flag & SOR_ZERO_INITIAL_GUESS) { 2012 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2013 for (i=0; i<m; i++) 2014 t[i] = b[i]; 2015 2016 for (i=0; i<m; i++){ 2017 d = *(aa + ai[i]); /* diag[i] */ 2018 v = aa + ai[i] + 1; 2019 vj = aj + ai[i] + 1; 2020 nz = ai[i+1] - ai[i] - 1; 2021 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2022 x[i] = omega*t[i]/d; 2023 while (nz--) t[*vj++] -= x[i]*(*v++); /* update rhs */ 2024 } 2025 } 2026 2027 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2028 if (!(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)){ 2029 t = b; 2030 } 2031 2032 for (i=m-1; i>=0; i--){ 2033 d = *(aa + ai[i]); 2034 v = aa + ai[i] + 1; 2035 vj = aj + ai[i] + 1; 2036 nz = ai[i+1] - ai[i] - 1; 2037 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2038 sum = t[i]; 2039 while (nz--) sum -= x[*vj++]*(*v++); 2040 x[i] = (1-omega)*x[i] + omega*sum/d; 2041 } 2042 t = a->relax_work; 2043 } 2044 its--; 2045 } 2046 2047 while (its--) { 2048 /* 2049 forward sweep: 2050 for i=0,...,m-1: 2051 sum[i] = (b[i] - U(i,:)x )/d[i]; 2052 x[i] = (1-omega)x[i] + omega*sum[i]; 2053 b = b - x[i]*U^T(i,:); 2054 2055 */ 2056 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2057 for (i=0; i<m; i++) 2058 t[i] = b[i]; 2059 2060 for (i=0; i<m; i++){ 2061 d = *(aa + ai[i]); /* diag[i] */ 2062 v = aa + ai[i] + 1; v1=v; 2063 vj = aj + ai[i] + 1; vj1=vj; 2064 nz = ai[i+1] - ai[i] - 1; nz1=nz; 2065 sum = t[i]; 2066 ierr = PetscLogFlops(4*nz-2);CHKERRQ(ierr); 2067 while (nz1--) sum -= (*v1++)*x[*vj1++]; 2068 x[i] = (1-omega)*x[i] + omega*sum/d; 2069 while (nz--) t[*vj++] -= x[i]*(*v++); 2070 } 2071 } 2072 2073 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2074 /* 2075 backward sweep: 2076 b = b - x[i]*U^T(i,:), i=0,...,n-2 2077 for i=m-1,...,0: 2078 sum[i] = (b[i] - U(i,:)x )/d[i]; 2079 x[i] = (1-omega)x[i] + omega*sum[i]; 2080 */ 2081 /* if there was a forward sweep done above then I thing the next two for loops are not needed */ 2082 for (i=0; i<m; i++) 2083 t[i] = b[i]; 2084 2085 for (i=0; i<m-1; i++){ /* update rhs */ 2086 v = aa + ai[i] + 1; 2087 vj = aj + ai[i] + 1; 2088 nz = ai[i+1] - ai[i] - 1; 2089 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2090 while (nz--) t[*vj++] -= x[i]*(*v++); 2091 } 2092 for (i=m-1; i>=0; i--){ 2093 d = *(aa + ai[i]); 2094 v = aa + ai[i] + 1; 2095 vj = aj + ai[i] + 1; 2096 nz = ai[i+1] - ai[i] - 1; 2097 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2098 sum = t[i]; 2099 while (nz--) sum -= x[*vj++]*(*v++); 2100 x[i] = (1-omega)*x[i] + omega*sum/d; 2101 } 2102 } 2103 } 2104 2105 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 2106 if (bb != xx) { 2107 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 2108 } 2109 PetscFunctionReturn(0); 2110 } 2111 2112 #undef __FUNCT__ 2113 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2114 /*@ 2115 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2116 (upper triangular entries in CSR format) provided by the user. 2117 2118 Collective on MPI_Comm 2119 2120 Input Parameters: 2121 + comm - must be an MPI communicator of size 1 2122 . bs - size of block 2123 . m - number of rows 2124 . n - number of columns 2125 . i - row indices 2126 . j - column indices 2127 - a - matrix values 2128 2129 Output Parameter: 2130 . mat - the matrix 2131 2132 Level: intermediate 2133 2134 Notes: 2135 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2136 once the matrix is destroyed 2137 2138 You cannot set new nonzero locations into this matrix, that will generate an error. 2139 2140 The i and j indices are 0 based 2141 2142 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2143 2144 @*/ 2145 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2146 { 2147 PetscErrorCode ierr; 2148 PetscInt ii; 2149 Mat_SeqSBAIJ *sbaij; 2150 2151 PetscFunctionBegin; 2152 if (bs != 1) SETERRQ1(PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2153 if (i[0]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2154 2155 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2156 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2157 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2158 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2159 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2160 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2161 2162 sbaij->i = i; 2163 sbaij->j = j; 2164 sbaij->a = a; 2165 sbaij->singlemalloc = PETSC_FALSE; 2166 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2167 sbaij->free_a = PETSC_FALSE; 2168 sbaij->free_ij = PETSC_FALSE; 2169 2170 for (ii=0; ii<m; ii++) { 2171 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2172 #if defined(PETSC_USE_DEBUG) 2173 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]); 2174 #endif 2175 } 2176 #if defined(PETSC_USE_DEBUG) 2177 for (ii=0; ii<sbaij->i[m]; ii++) { 2178 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2179 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2180 } 2181 #endif 2182 2183 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2184 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2185 PetscFunctionReturn(0); 2186 } 2187 2188 2189 2190 2191 2192