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