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