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