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 #undef __FUNCT__ 941 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 942 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,MatFactorInfo *info) 943 { 944 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 945 Mat outA; 946 PetscErrorCode ierr; 947 PetscTruth row_identity; 948 949 PetscFunctionBegin; 950 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 951 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 952 if (!row_identity) SETERRQ(PETSC_ERR_SUP,"Matrix reordering is not supported"); 953 if (inA->rmap.bs != 1) SETERRQ1(PETSC_ERR_SUP,"Matrix block size %D is not supported",inA->rmap.bs); /* Need to replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR()! */ 954 955 outA = inA; 956 inA->factor = FACTOR_CHOLESKY; 957 958 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 959 /* 960 Blocksize < 8 have a special faster factorization/solver 961 for ICC(0) factorization with natural ordering 962 */ 963 switch (inA->rmap.bs){ /* Note: row_identity = PETSC_TRUE! */ 964 case 1: 965 inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering; 966 inA->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering; 967 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering; 968 inA->ops->solves = MatSolves_SeqSBAIJ_1; 969 ierr = PetscInfo(inA,"Using special in-place natural ordering solvetrans BS=1\n");CHKERRQ(ierr); 970 break; 971 case 2: 972 inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 973 inA->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering; 974 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering; 975 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=2\n");CHKERRQ(ierr); 976 break; 977 case 3: 978 inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 979 inA->ops->solve = MatSolve_SeqSBAIJ_3_NaturalOrdering; 980 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_3_NaturalOrdering; 981 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=3\n");CHKERRQ(ierr); 982 break; 983 case 4: 984 inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 985 inA->ops->solve = MatSolve_SeqSBAIJ_4_NaturalOrdering; 986 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_4_NaturalOrdering; 987 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=4\n");CHKERRQ(ierr); 988 break; 989 case 5: 990 inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 991 inA->ops->solve = MatSolve_SeqSBAIJ_5_NaturalOrdering; 992 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_5_NaturalOrdering; 993 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=5\n");CHKERRQ(ierr); 994 break; 995 case 6: 996 inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 997 inA->ops->solve = MatSolve_SeqSBAIJ_6_NaturalOrdering; 998 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_6_NaturalOrdering; 999 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=6\n");CHKERRQ(ierr); 1000 break; 1001 case 7: 1002 inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1003 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_7_NaturalOrdering; 1004 inA->ops->solve = MatSolve_SeqSBAIJ_7_NaturalOrdering; 1005 ierr = PetscInfo(inA,"Using special in-place natural ordering factor and solve BS=7\n");CHKERRQ(ierr); 1006 break; 1007 default: 1008 inA->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1009 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_N_NaturalOrdering; 1010 inA->ops->solve = MatSolve_SeqSBAIJ_N_NaturalOrdering; 1011 break; 1012 } 1013 1014 a->row = row; 1015 a->col = row; 1016 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1017 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1018 1019 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1020 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1021 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1022 1023 if (!a->solve_work) { 1024 ierr = PetscMalloc((inA->rmap.N+inA->rmap.bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1025 ierr = PetscLogObjectMemory(inA,(inA->rmap.N+inA->rmap.bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1026 } 1027 1028 ierr = MatCholeskyFactorNumeric(inA,info,&outA);CHKERRQ(ierr); 1029 PetscFunctionReturn(0); 1030 } 1031 1032 EXTERN_C_BEGIN 1033 #undef __FUNCT__ 1034 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1035 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1036 { 1037 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1038 PetscInt i,nz,n; 1039 1040 PetscFunctionBegin; 1041 nz = baij->maxnz; 1042 n = mat->cmap.n; 1043 for (i=0; i<nz; i++) { 1044 baij->j[i] = indices[i]; 1045 } 1046 baij->nz = nz; 1047 for (i=0; i<n; i++) { 1048 baij->ilen[i] = baij->imax[i]; 1049 } 1050 PetscFunctionReturn(0); 1051 } 1052 EXTERN_C_END 1053 1054 #undef __FUNCT__ 1055 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1056 /*@ 1057 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1058 in the matrix. 1059 1060 Input Parameters: 1061 + mat - the SeqSBAIJ matrix 1062 - indices - the column indices 1063 1064 Level: advanced 1065 1066 Notes: 1067 This can be called if you have precomputed the nonzero structure of the 1068 matrix and want to provide it to the matrix object to improve the performance 1069 of the MatSetValues() operation. 1070 1071 You MUST have set the correct numbers of nonzeros per row in the call to 1072 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1073 1074 MUST be called before any calls to MatSetValues() 1075 1076 .seealso: MatCreateSeqSBAIJ 1077 @*/ 1078 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1079 { 1080 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 1081 1082 PetscFunctionBegin; 1083 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 1084 PetscValidPointer(indices,2); 1085 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 1086 if (f) { 1087 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1088 } else { 1089 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 1090 } 1091 PetscFunctionReturn(0); 1092 } 1093 1094 #undef __FUNCT__ 1095 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1096 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1097 { 1098 PetscErrorCode ierr; 1099 1100 PetscFunctionBegin; 1101 /* If the two matrices have the same copy implementation, use fast copy. */ 1102 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1103 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1104 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1105 1106 if (a->i[A->rmap.N] != b->i[B->rmap.N]) { 1107 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1108 } 1109 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap.N])*sizeof(PetscScalar));CHKERRQ(ierr); 1110 } else { 1111 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1112 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1113 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1114 } 1115 PetscFunctionReturn(0); 1116 } 1117 1118 #undef __FUNCT__ 1119 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1120 PetscErrorCode MatSetUpPreallocation_SeqSBAIJ(Mat A) 1121 { 1122 PetscErrorCode ierr; 1123 1124 PetscFunctionBegin; 1125 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,PetscMax(A->rmap.bs,1),PETSC_DEFAULT,0);CHKERRQ(ierr); 1126 PetscFunctionReturn(0); 1127 } 1128 1129 #undef __FUNCT__ 1130 #define __FUNCT__ "MatGetArray_SeqSBAIJ" 1131 PetscErrorCode MatGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1132 { 1133 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1134 PetscFunctionBegin; 1135 *array = a->a; 1136 PetscFunctionReturn(0); 1137 } 1138 1139 #undef __FUNCT__ 1140 #define __FUNCT__ "MatRestoreArray_SeqSBAIJ" 1141 PetscErrorCode MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1142 { 1143 PetscFunctionBegin; 1144 PetscFunctionReturn(0); 1145 } 1146 1147 #include "petscblaslapack.h" 1148 #undef __FUNCT__ 1149 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1150 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1151 { 1152 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ *)X->data, *y=(Mat_SeqSBAIJ *)Y->data; 1153 PetscErrorCode ierr; 1154 PetscInt i,bs=Y->rmap.bs,bs2,j; 1155 PetscBLASInt bnz = (PetscBLASInt)x->nz,one = 1; 1156 1157 PetscFunctionBegin; 1158 if (str == SAME_NONZERO_PATTERN) { 1159 PetscScalar alpha = a; 1160 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 1161 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1162 if (y->xtoy && y->XtoY != X) { 1163 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 1164 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 1165 } 1166 if (!y->xtoy) { /* get xtoy */ 1167 ierr = MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 1168 y->XtoY = X; 1169 } 1170 bs2 = bs*bs; 1171 for (i=0; i<x->nz; i++) { 1172 j = 0; 1173 while (j < bs2){ 1174 y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]); 1175 j++; 1176 } 1177 } 1178 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); 1179 } else { 1180 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1181 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1182 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1183 } 1184 PetscFunctionReturn(0); 1185 } 1186 1187 #undef __FUNCT__ 1188 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1189 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg) 1190 { 1191 PetscFunctionBegin; 1192 *flg = PETSC_TRUE; 1193 PetscFunctionReturn(0); 1194 } 1195 1196 #undef __FUNCT__ 1197 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1198 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscTruth *flg) 1199 { 1200 PetscFunctionBegin; 1201 *flg = PETSC_TRUE; 1202 PetscFunctionReturn(0); 1203 } 1204 1205 #undef __FUNCT__ 1206 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1207 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscTruth *flg) 1208 { 1209 PetscFunctionBegin; 1210 *flg = PETSC_FALSE; 1211 PetscFunctionReturn(0); 1212 } 1213 1214 #undef __FUNCT__ 1215 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1216 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1217 { 1218 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1219 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1220 PetscScalar *aa = a->a; 1221 1222 PetscFunctionBegin; 1223 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1224 PetscFunctionReturn(0); 1225 } 1226 1227 #undef __FUNCT__ 1228 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1229 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1230 { 1231 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1232 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1233 PetscScalar *aa = a->a; 1234 1235 PetscFunctionBegin; 1236 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1237 PetscFunctionReturn(0); 1238 } 1239 1240 /* -------------------------------------------------------------------*/ 1241 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1242 MatGetRow_SeqSBAIJ, 1243 MatRestoreRow_SeqSBAIJ, 1244 MatMult_SeqSBAIJ_N, 1245 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1246 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1247 MatMultAdd_SeqSBAIJ_N, 1248 MatSolve_SeqSBAIJ_N, 1249 0, 1250 0, 1251 /*10*/ 0, 1252 0, 1253 MatCholeskyFactor_SeqSBAIJ, 1254 MatRelax_SeqSBAIJ, 1255 MatTranspose_SeqSBAIJ, 1256 /*15*/ MatGetInfo_SeqSBAIJ, 1257 MatEqual_SeqSBAIJ, 1258 MatGetDiagonal_SeqSBAIJ, 1259 MatDiagonalScale_SeqSBAIJ, 1260 MatNorm_SeqSBAIJ, 1261 /*20*/ 0, 1262 MatAssemblyEnd_SeqSBAIJ, 1263 0, 1264 MatSetOption_SeqSBAIJ, 1265 MatZeroEntries_SeqSBAIJ, 1266 /*25*/ 0, 1267 0, 1268 0, 1269 MatCholeskyFactorSymbolic_SeqSBAIJ, 1270 MatCholeskyFactorNumeric_SeqSBAIJ_N, 1271 /*30*/ MatSetUpPreallocation_SeqSBAIJ, 1272 0, 1273 MatICCFactorSymbolic_SeqSBAIJ, 1274 MatGetArray_SeqSBAIJ, 1275 MatRestoreArray_SeqSBAIJ, 1276 /*35*/ MatDuplicate_SeqSBAIJ, 1277 0, 1278 0, 1279 0, 1280 MatICCFactor_SeqSBAIJ, 1281 /*40*/ MatAXPY_SeqSBAIJ, 1282 MatGetSubMatrices_SeqSBAIJ, 1283 MatIncreaseOverlap_SeqSBAIJ, 1284 MatGetValues_SeqSBAIJ, 1285 MatCopy_SeqSBAIJ, 1286 /*45*/ 0, 1287 MatScale_SeqSBAIJ, 1288 0, 1289 0, 1290 0, 1291 /*50*/ 0, 1292 MatGetRowIJ_SeqSBAIJ, 1293 MatRestoreRowIJ_SeqSBAIJ, 1294 0, 1295 0, 1296 /*55*/ 0, 1297 0, 1298 0, 1299 0, 1300 MatSetValuesBlocked_SeqSBAIJ, 1301 /*60*/ MatGetSubMatrix_SeqSBAIJ, 1302 0, 1303 0, 1304 0, 1305 0, 1306 /*65*/ 0, 1307 0, 1308 0, 1309 0, 1310 0, 1311 /*70*/ MatGetRowMax_SeqSBAIJ, 1312 0, 1313 0, 1314 0, 1315 0, 1316 /*75*/ 0, 1317 0, 1318 0, 1319 0, 1320 0, 1321 /*80*/ 0, 1322 0, 1323 0, 1324 #if !defined(PETSC_USE_COMPLEX) 1325 MatGetInertia_SeqSBAIJ, 1326 #else 1327 0, 1328 #endif 1329 MatLoad_SeqSBAIJ, 1330 /*85*/ MatIsSymmetric_SeqSBAIJ, 1331 MatIsHermitian_SeqSBAIJ, 1332 MatIsStructurallySymmetric_SeqSBAIJ, 1333 0, 1334 0, 1335 /*90*/ 0, 1336 0, 1337 0, 1338 0, 1339 0, 1340 /*95*/ 0, 1341 0, 1342 0, 1343 0, 1344 0, 1345 /*100*/0, 1346 0, 1347 0, 1348 0, 1349 0, 1350 /*105*/0, 1351 MatRealPart_SeqSBAIJ, 1352 MatImaginaryPart_SeqSBAIJ, 1353 MatGetRowUpperTriangular_SeqSBAIJ, 1354 MatRestoreRowUpperTriangular_SeqSBAIJ 1355 }; 1356 1357 EXTERN_C_BEGIN 1358 #undef __FUNCT__ 1359 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1360 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqSBAIJ(Mat mat) 1361 { 1362 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1363 PetscInt nz = aij->i[mat->rmap.N]*mat->rmap.bs*aij->bs2; 1364 PetscErrorCode ierr; 1365 1366 PetscFunctionBegin; 1367 if (aij->nonew != 1) { 1368 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1369 } 1370 1371 /* allocate space for values if not already there */ 1372 if (!aij->saved_values) { 1373 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1374 } 1375 1376 /* copy values over */ 1377 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1378 PetscFunctionReturn(0); 1379 } 1380 EXTERN_C_END 1381 1382 EXTERN_C_BEGIN 1383 #undef __FUNCT__ 1384 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1385 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqSBAIJ(Mat mat) 1386 { 1387 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1388 PetscErrorCode ierr; 1389 PetscInt nz = aij->i[mat->rmap.N]*mat->rmap.bs*aij->bs2; 1390 1391 PetscFunctionBegin; 1392 if (aij->nonew != 1) { 1393 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1394 } 1395 if (!aij->saved_values) { 1396 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1397 } 1398 1399 /* copy values over */ 1400 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1401 PetscFunctionReturn(0); 1402 } 1403 EXTERN_C_END 1404 1405 EXTERN_C_BEGIN 1406 #undef __FUNCT__ 1407 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1408 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1409 { 1410 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1411 PetscErrorCode ierr; 1412 PetscInt i,mbs,bs2; 1413 PetscTruth skipallocation = PETSC_FALSE,flg; 1414 1415 PetscFunctionBegin; 1416 B->preallocated = PETSC_TRUE; 1417 ierr = PetscOptionsGetInt(B->prefix,"-mat_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1418 B->rmap.bs = B->cmap.bs = bs; 1419 ierr = PetscMapInitialize(B->comm,&B->rmap);CHKERRQ(ierr); 1420 ierr = PetscMapInitialize(B->comm,&B->cmap);CHKERRQ(ierr); 1421 1422 mbs = B->rmap.N/bs; 1423 bs2 = bs*bs; 1424 1425 if (mbs*bs != B->rmap.N) { 1426 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1427 } 1428 1429 if (nz == MAT_SKIP_ALLOCATION) { 1430 skipallocation = PETSC_TRUE; 1431 nz = 0; 1432 } 1433 1434 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1435 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1436 if (nnz) { 1437 for (i=0; i<mbs; i++) { 1438 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]); 1439 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); 1440 } 1441 } 1442 1443 ierr = PetscOptionsHasName(B->prefix,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1444 if (!flg) { 1445 switch (bs) { 1446 case 1: 1447 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1448 B->ops->solve = MatSolve_SeqSBAIJ_1; 1449 B->ops->solves = MatSolves_SeqSBAIJ_1; 1450 B->ops->solvetranspose = MatSolve_SeqSBAIJ_1; 1451 B->ops->mult = MatMult_SeqSBAIJ_1; 1452 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1453 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1454 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1455 break; 1456 case 2: 1457 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1458 B->ops->solve = MatSolve_SeqSBAIJ_2; 1459 B->ops->solvetranspose = MatSolve_SeqSBAIJ_2; 1460 B->ops->mult = MatMult_SeqSBAIJ_2; 1461 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1462 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1463 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1464 break; 1465 case 3: 1466 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1467 B->ops->solve = MatSolve_SeqSBAIJ_3; 1468 B->ops->solvetranspose = MatSolve_SeqSBAIJ_3; 1469 B->ops->mult = MatMult_SeqSBAIJ_3; 1470 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1471 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1472 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1473 break; 1474 case 4: 1475 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1476 B->ops->solve = MatSolve_SeqSBAIJ_4; 1477 B->ops->solvetranspose = MatSolve_SeqSBAIJ_4; 1478 B->ops->mult = MatMult_SeqSBAIJ_4; 1479 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1480 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1481 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1482 break; 1483 case 5: 1484 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1485 B->ops->solve = MatSolve_SeqSBAIJ_5; 1486 B->ops->solvetranspose = MatSolve_SeqSBAIJ_5; 1487 B->ops->mult = MatMult_SeqSBAIJ_5; 1488 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1489 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1490 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1491 break; 1492 case 6: 1493 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1494 B->ops->solve = MatSolve_SeqSBAIJ_6; 1495 B->ops->solvetranspose = MatSolve_SeqSBAIJ_6; 1496 B->ops->mult = MatMult_SeqSBAIJ_6; 1497 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1498 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1499 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1500 break; 1501 case 7: 1502 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1503 B->ops->solve = MatSolve_SeqSBAIJ_7; 1504 B->ops->solvetranspose = MatSolve_SeqSBAIJ_7; 1505 B->ops->mult = MatMult_SeqSBAIJ_7; 1506 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1507 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1508 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1509 break; 1510 } 1511 } 1512 1513 b->mbs = mbs; 1514 b->nbs = mbs; 1515 if (!skipallocation) { 1516 /* b->ilen will count nonzeros in each block row so far. */ 1517 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1518 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1519 if (!nnz) { 1520 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1521 else if (nz <= 0) nz = 1; 1522 for (i=0; i<mbs; i++) { 1523 b->imax[i] = nz; 1524 } 1525 nz = nz*mbs; /* total nz */ 1526 } else { 1527 nz = 0; 1528 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1529 } 1530 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1531 1532 /* allocate the matrix space */ 1533 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap.N+1,PetscInt,&b->i);CHKERRQ(ierr); 1534 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1535 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1536 b->singlemalloc = PETSC_TRUE; 1537 1538 /* pointer to beginning of each row */ 1539 b->i[0] = 0; 1540 for (i=1; i<mbs+1; i++) { 1541 b->i[i] = b->i[i-1] + b->imax[i-1]; 1542 } 1543 b->free_a = PETSC_TRUE; 1544 b->free_ij = PETSC_TRUE; 1545 } else { 1546 b->free_a = PETSC_FALSE; 1547 b->free_ij = PETSC_FALSE; 1548 } 1549 1550 B->rmap.bs = bs; 1551 b->bs2 = bs2; 1552 b->nz = 0; 1553 b->maxnz = nz*bs2; 1554 1555 b->inew = 0; 1556 b->jnew = 0; 1557 b->anew = 0; 1558 b->a2anew = 0; 1559 b->permute = PETSC_FALSE; 1560 PetscFunctionReturn(0); 1561 } 1562 EXTERN_C_END 1563 1564 EXTERN_C_BEGIN 1565 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1566 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1567 EXTERN_C_END 1568 1569 /*MC 1570 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1571 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1572 1573 Options Database Keys: 1574 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1575 1576 Level: beginner 1577 1578 .seealso: MatCreateSeqSBAIJ 1579 M*/ 1580 1581 EXTERN_C_BEGIN 1582 #undef __FUNCT__ 1583 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1584 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqSBAIJ(Mat B) 1585 { 1586 Mat_SeqSBAIJ *b; 1587 PetscErrorCode ierr; 1588 PetscMPIInt size; 1589 PetscTruth flg; 1590 1591 PetscFunctionBegin; 1592 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1593 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1594 1595 ierr = PetscNew(Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1596 B->data = (void*)b; 1597 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1598 B->ops->destroy = MatDestroy_SeqSBAIJ; 1599 B->ops->view = MatView_SeqSBAIJ; 1600 B->factor = 0; 1601 B->mapping = 0; 1602 b->row = 0; 1603 b->icol = 0; 1604 b->reallocs = 0; 1605 b->saved_values = 0; 1606 1607 1608 b->sorted = PETSC_FALSE; 1609 b->roworiented = PETSC_TRUE; 1610 b->nonew = 0; 1611 b->diag = 0; 1612 b->solve_work = 0; 1613 b->mult_work = 0; 1614 B->spptr = 0; 1615 b->keepzeroedrows = PETSC_FALSE; 1616 b->xtoy = 0; 1617 b->XtoY = 0; 1618 1619 b->inew = 0; 1620 b->jnew = 0; 1621 b->anew = 0; 1622 b->a2anew = 0; 1623 b->permute = PETSC_FALSE; 1624 1625 b->ignore_ltriangular = PETSC_FALSE; 1626 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_ignore_lower_triangular",&flg);CHKERRQ(ierr); 1627 if (flg) b->ignore_ltriangular = PETSC_TRUE; 1628 1629 b->getrow_utriangular = PETSC_FALSE; 1630 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_getrow_uppertriangular",&flg);CHKERRQ(ierr); 1631 if (flg) b->getrow_utriangular = PETSC_TRUE; 1632 1633 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1634 "MatStoreValues_SeqSBAIJ", 1635 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1636 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1637 "MatRetrieveValues_SeqSBAIJ", 1638 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1639 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1640 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1641 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1642 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1643 "MatConvert_SeqSBAIJ_SeqAIJ", 1644 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1645 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1646 "MatConvert_SeqSBAIJ_SeqBAIJ", 1647 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1648 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1649 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1650 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1651 1652 B->symmetric = PETSC_TRUE; 1653 B->structurally_symmetric = PETSC_TRUE; 1654 B->symmetric_set = PETSC_TRUE; 1655 B->structurally_symmetric_set = PETSC_TRUE; 1656 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1657 PetscFunctionReturn(0); 1658 } 1659 EXTERN_C_END 1660 1661 #undef __FUNCT__ 1662 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1663 /*@C 1664 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1665 compressed row) format. For good matrix assembly performance the 1666 user should preallocate the matrix storage by setting the parameter nz 1667 (or the array nnz). By setting these parameters accurately, performance 1668 during matrix assembly can be increased by more than a factor of 50. 1669 1670 Collective on Mat 1671 1672 Input Parameters: 1673 + A - the symmetric matrix 1674 . bs - size of block 1675 . nz - number of block nonzeros per block row (same for all rows) 1676 - nnz - array containing the number of block nonzeros in the upper triangular plus 1677 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1678 1679 Options Database Keys: 1680 . -mat_no_unroll - uses code that does not unroll the loops in the 1681 block calculations (much slower) 1682 . -mat_block_size - size of the blocks to use 1683 1684 Level: intermediate 1685 1686 Notes: 1687 Specify the preallocated storage with either nz or nnz (not both). 1688 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1689 allocation. For additional details, see the users manual chapter on 1690 matrices. 1691 1692 If the nnz parameter is given then the nz parameter is ignored 1693 1694 1695 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1696 @*/ 1697 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1698 { 1699 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]); 1700 1701 PetscFunctionBegin; 1702 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1703 if (f) { 1704 ierr = (*f)(B,bs,nz,nnz);CHKERRQ(ierr); 1705 } 1706 PetscFunctionReturn(0); 1707 } 1708 1709 #undef __FUNCT__ 1710 #define __FUNCT__ "MatCreateSeqSBAIJ" 1711 /*@C 1712 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1713 compressed row) format. For good matrix assembly performance the 1714 user should preallocate the matrix storage by setting the parameter nz 1715 (or the array nnz). By setting these parameters accurately, performance 1716 during matrix assembly can be increased by more than a factor of 50. 1717 1718 Collective on MPI_Comm 1719 1720 Input Parameters: 1721 + comm - MPI communicator, set to PETSC_COMM_SELF 1722 . bs - size of block 1723 . m - number of rows, or number of columns 1724 . nz - number of block nonzeros per block row (same for all rows) 1725 - nnz - array containing the number of block nonzeros in the upper triangular plus 1726 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1727 1728 Output Parameter: 1729 . A - the symmetric matrix 1730 1731 Options Database Keys: 1732 . -mat_no_unroll - uses code that does not unroll the loops in the 1733 block calculations (much slower) 1734 . -mat_block_size - size of the blocks to use 1735 1736 Level: intermediate 1737 1738 Notes: 1739 1740 Specify the preallocated storage with either nz or nnz (not both). 1741 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1742 allocation. For additional details, see the users manual chapter on 1743 matrices. 1744 1745 If the nnz parameter is given then the nz parameter is ignored 1746 1747 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1748 @*/ 1749 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 1750 { 1751 PetscErrorCode ierr; 1752 1753 PetscFunctionBegin; 1754 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1755 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1756 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1757 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 1758 PetscFunctionReturn(0); 1759 } 1760 1761 #undef __FUNCT__ 1762 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1763 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1764 { 1765 Mat C; 1766 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1767 PetscErrorCode ierr; 1768 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1769 1770 PetscFunctionBegin; 1771 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1772 1773 *B = 0; 1774 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1775 ierr = MatSetSizes(C,A->rmap.N,A->cmap.n,A->rmap.N,A->cmap.n);CHKERRQ(ierr); 1776 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1777 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1778 c = (Mat_SeqSBAIJ*)C->data; 1779 1780 C->preallocated = PETSC_TRUE; 1781 C->factor = A->factor; 1782 c->row = 0; 1783 c->icol = 0; 1784 c->saved_values = 0; 1785 c->keepzeroedrows = a->keepzeroedrows; 1786 C->assembled = PETSC_TRUE; 1787 1788 ierr = PetscMapCopy(A->comm,&A->rmap,&C->rmap);CHKERRQ(ierr); 1789 ierr = PetscMapCopy(A->comm,&A->cmap,&C->cmap);CHKERRQ(ierr); 1790 c->bs2 = a->bs2; 1791 c->mbs = a->mbs; 1792 c->nbs = a->nbs; 1793 1794 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 1795 for (i=0; i<mbs; i++) { 1796 c->imax[i] = a->imax[i]; 1797 c->ilen[i] = a->ilen[i]; 1798 } 1799 1800 /* allocate the matrix space */ 1801 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 1802 c->singlemalloc = PETSC_TRUE; 1803 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1804 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 1805 if (mbs > 0) { 1806 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1807 if (cpvalues == MAT_COPY_VALUES) { 1808 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1809 } else { 1810 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1811 } 1812 } 1813 1814 c->sorted = a->sorted; 1815 c->roworiented = a->roworiented; 1816 c->nonew = a->nonew; 1817 1818 if (a->diag) { 1819 ierr = PetscMalloc((mbs+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 1820 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1821 for (i=0; i<mbs; i++) { 1822 c->diag[i] = a->diag[i]; 1823 } 1824 } else c->diag = 0; 1825 c->nz = a->nz; 1826 c->maxnz = a->maxnz; 1827 c->solve_work = 0; 1828 c->mult_work = 0; 1829 c->free_a = PETSC_TRUE; 1830 c->free_ij = PETSC_TRUE; 1831 *B = C; 1832 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1833 PetscFunctionReturn(0); 1834 } 1835 1836 #undef __FUNCT__ 1837 #define __FUNCT__ "MatLoad_SeqSBAIJ" 1838 PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, MatType type,Mat *A) 1839 { 1840 Mat_SeqSBAIJ *a; 1841 Mat B; 1842 PetscErrorCode ierr; 1843 int fd; 1844 PetscMPIInt size; 1845 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 1846 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 1847 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows; 1848 PetscInt *masked,nmask,tmp,bs2,ishift; 1849 PetscScalar *aa; 1850 MPI_Comm comm = ((PetscObject)viewer)->comm; 1851 1852 PetscFunctionBegin; 1853 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1854 bs2 = bs*bs; 1855 1856 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1857 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1858 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1859 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1860 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1861 M = header[1]; N = header[2]; nz = header[3]; 1862 1863 if (header[3] < 0) { 1864 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 1865 } 1866 1867 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1868 1869 /* 1870 This code adds extra rows to make sure the number of rows is 1871 divisible by the blocksize 1872 */ 1873 mbs = M/bs; 1874 extra_rows = bs - M + bs*(mbs); 1875 if (extra_rows == bs) extra_rows = 0; 1876 else mbs++; 1877 if (extra_rows) { 1878 ierr = PetscInfo(0,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 1879 } 1880 1881 /* read in row lengths */ 1882 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 1883 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1884 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1885 1886 /* read in column indices */ 1887 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 1888 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1889 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1890 1891 /* loop over row lengths determining block row lengths */ 1892 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 1893 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 1894 ierr = PetscMalloc(2*mbs*sizeof(PetscInt),&mask);CHKERRQ(ierr); 1895 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 1896 masked = mask + mbs; 1897 rowcount = 0; nzcount = 0; 1898 for (i=0; i<mbs; i++) { 1899 nmask = 0; 1900 for (j=0; j<bs; j++) { 1901 kmax = rowlengths[rowcount]; 1902 for (k=0; k<kmax; k++) { 1903 tmp = jj[nzcount++]/bs; /* block col. index */ 1904 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 1905 } 1906 rowcount++; 1907 } 1908 s_browlengths[i] += nmask; 1909 1910 /* zero out the mask elements we set */ 1911 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1912 } 1913 1914 /* create our matrix */ 1915 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 1916 ierr = MatSetSizes(B,M+extra_rows,N+extra_rows,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 1917 ierr = MatSetType(B,type);CHKERRQ(ierr); 1918 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,s_browlengths);CHKERRQ(ierr); 1919 a = (Mat_SeqSBAIJ*)B->data; 1920 1921 /* set matrix "i" values */ 1922 a->i[0] = 0; 1923 for (i=1; i<= mbs; i++) { 1924 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 1925 a->ilen[i-1] = s_browlengths[i-1]; 1926 } 1927 a->nz = a->i[mbs]; 1928 1929 /* read in nonzero values */ 1930 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 1931 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1932 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1933 1934 /* set "a" and "j" values into matrix */ 1935 nzcount = 0; jcount = 0; 1936 for (i=0; i<mbs; i++) { 1937 nzcountb = nzcount; 1938 nmask = 0; 1939 for (j=0; j<bs; j++) { 1940 kmax = rowlengths[i*bs+j]; 1941 for (k=0; k<kmax; k++) { 1942 tmp = jj[nzcount++]/bs; /* block col. index */ 1943 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 1944 } 1945 } 1946 /* sort the masked values */ 1947 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1948 1949 /* set "j" values into matrix */ 1950 maskcount = 1; 1951 for (j=0; j<nmask; j++) { 1952 a->j[jcount++] = masked[j]; 1953 mask[masked[j]] = maskcount++; 1954 } 1955 1956 /* set "a" values into matrix */ 1957 ishift = bs2*a->i[i]; 1958 for (j=0; j<bs; j++) { 1959 kmax = rowlengths[i*bs+j]; 1960 for (k=0; k<kmax; k++) { 1961 tmp = jj[nzcountb]/bs ; /* block col. index */ 1962 if (tmp >= i){ 1963 block = mask[tmp] - 1; 1964 point = jj[nzcountb] - bs*tmp; 1965 idx = ishift + bs2*block + j + bs*point; 1966 a->a[idx] = aa[nzcountb]; 1967 } 1968 nzcountb++; 1969 } 1970 } 1971 /* zero out the mask elements we set */ 1972 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1973 } 1974 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 1975 1976 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1977 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 1978 ierr = PetscFree(aa);CHKERRQ(ierr); 1979 ierr = PetscFree(jj);CHKERRQ(ierr); 1980 ierr = PetscFree(mask);CHKERRQ(ierr); 1981 1982 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1983 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1984 ierr = MatView_Private(B);CHKERRQ(ierr); 1985 *A = B; 1986 PetscFunctionReturn(0); 1987 } 1988 1989 #undef __FUNCT__ 1990 #define __FUNCT__ "MatRelax_SeqSBAIJ" 1991 PetscErrorCode MatRelax_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1992 { 1993 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1994 MatScalar *aa=a->a,*v,*v1; 1995 PetscScalar *x,*b,*t,sum,d; 1996 PetscErrorCode ierr; 1997 PetscInt m=a->mbs,bs=A->rmap.bs,*ai=a->i,*aj=a->j; 1998 PetscInt nz,nz1,*vj,*vj1,i; 1999 2000 PetscFunctionBegin; 2001 if (flag & SOR_EISENSTAT) SETERRQ(PETSC_ERR_SUP,"No support yet for Eisenstat"); 2002 2003 its = its*lits; 2004 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 2005 2006 if (bs > 1) 2007 SETERRQ(PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented"); 2008 2009 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 2010 if (xx != bb) { 2011 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 2012 } else { 2013 b = x; 2014 } 2015 2016 if (!a->relax_work) { 2017 ierr = PetscMalloc(m*sizeof(PetscScalar),&a->relax_work);CHKERRQ(ierr); 2018 } 2019 t = a->relax_work; 2020 2021 if (flag & SOR_ZERO_INITIAL_GUESS) { 2022 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2023 for (i=0; i<m; i++) 2024 t[i] = b[i]; 2025 2026 for (i=0; i<m; i++){ 2027 d = *(aa + ai[i]); /* diag[i] */ 2028 v = aa + ai[i] + 1; 2029 vj = aj + ai[i] + 1; 2030 nz = ai[i+1] - ai[i] - 1; 2031 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2032 x[i] = omega*t[i]/d; 2033 while (nz--) t[*vj++] -= x[i]*(*v++); /* update rhs */ 2034 } 2035 } 2036 2037 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2038 if (!(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)){ 2039 t = b; 2040 } 2041 2042 for (i=m-1; i>=0; i--){ 2043 d = *(aa + ai[i]); 2044 v = aa + ai[i] + 1; 2045 vj = aj + ai[i] + 1; 2046 nz = ai[i+1] - ai[i] - 1; 2047 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2048 sum = t[i]; 2049 while (nz--) sum -= x[*vj++]*(*v++); 2050 x[i] = (1-omega)*x[i] + omega*sum/d; 2051 } 2052 t = a->relax_work; 2053 } 2054 its--; 2055 } 2056 2057 while (its--) { 2058 /* 2059 forward sweep: 2060 for i=0,...,m-1: 2061 sum[i] = (b[i] - U(i,:)x )/d[i]; 2062 x[i] = (1-omega)x[i] + omega*sum[i]; 2063 b = b - x[i]*U^T(i,:); 2064 2065 */ 2066 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2067 for (i=0; i<m; i++) 2068 t[i] = b[i]; 2069 2070 for (i=0; i<m; i++){ 2071 d = *(aa + ai[i]); /* diag[i] */ 2072 v = aa + ai[i] + 1; v1=v; 2073 vj = aj + ai[i] + 1; vj1=vj; 2074 nz = ai[i+1] - ai[i] - 1; nz1=nz; 2075 sum = t[i]; 2076 ierr = PetscLogFlops(4*nz-2);CHKERRQ(ierr); 2077 while (nz1--) sum -= (*v1++)*x[*vj1++]; 2078 x[i] = (1-omega)*x[i] + omega*sum/d; 2079 while (nz--) t[*vj++] -= x[i]*(*v++); 2080 } 2081 } 2082 2083 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2084 /* 2085 backward sweep: 2086 b = b - x[i]*U^T(i,:), i=0,...,n-2 2087 for i=m-1,...,0: 2088 sum[i] = (b[i] - U(i,:)x )/d[i]; 2089 x[i] = (1-omega)x[i] + omega*sum[i]; 2090 */ 2091 /* if there was a forward sweep done above then I thing the next two for loops are not needed */ 2092 for (i=0; i<m; i++) 2093 t[i] = b[i]; 2094 2095 for (i=0; i<m-1; i++){ /* update rhs */ 2096 v = aa + ai[i] + 1; 2097 vj = aj + ai[i] + 1; 2098 nz = ai[i+1] - ai[i] - 1; 2099 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2100 while (nz--) t[*vj++] -= x[i]*(*v++); 2101 } 2102 for (i=m-1; i>=0; i--){ 2103 d = *(aa + ai[i]); 2104 v = aa + ai[i] + 1; 2105 vj = aj + ai[i] + 1; 2106 nz = ai[i+1] - ai[i] - 1; 2107 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2108 sum = t[i]; 2109 while (nz--) sum -= x[*vj++]*(*v++); 2110 x[i] = (1-omega)*x[i] + omega*sum/d; 2111 } 2112 } 2113 } 2114 2115 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 2116 if (bb != xx) { 2117 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 2118 } 2119 PetscFunctionReturn(0); 2120 } 2121 2122 #undef __FUNCT__ 2123 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2124 /*@ 2125 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2126 (upper triangular entries in CSR format) provided by the user. 2127 2128 Collective on MPI_Comm 2129 2130 Input Parameters: 2131 + comm - must be an MPI communicator of size 1 2132 . bs - size of block 2133 . m - number of rows 2134 . n - number of columns 2135 . i - row indices 2136 . j - column indices 2137 - a - matrix values 2138 2139 Output Parameter: 2140 . mat - the matrix 2141 2142 Level: intermediate 2143 2144 Notes: 2145 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2146 once the matrix is destroyed 2147 2148 You cannot set new nonzero locations into this matrix, that will generate an error. 2149 2150 The i and j indices are 0 based 2151 2152 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2153 2154 @*/ 2155 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2156 { 2157 PetscErrorCode ierr; 2158 PetscInt ii; 2159 Mat_SeqSBAIJ *sbaij; 2160 2161 PetscFunctionBegin; 2162 if (bs != 1) SETERRQ1(PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2163 if (i[0]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2164 2165 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2166 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2167 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2168 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2169 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2170 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2171 2172 sbaij->i = i; 2173 sbaij->j = j; 2174 sbaij->a = a; 2175 sbaij->singlemalloc = PETSC_FALSE; 2176 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2177 sbaij->free_a = PETSC_FALSE; 2178 sbaij->free_ij = PETSC_FALSE; 2179 2180 for (ii=0; ii<m; ii++) { 2181 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2182 #if defined(PETSC_USE_DEBUG) 2183 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]); 2184 #endif 2185 } 2186 #if defined(PETSC_USE_DEBUG) 2187 for (ii=0; ii<sbaij->i[m]; ii++) { 2188 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2189 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2190 } 2191 #endif 2192 2193 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2194 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2195 PetscFunctionReturn(0); 2196 } 2197 2198 2199 2200 2201 2202