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