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