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