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