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