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) stepval = (n-1)*bs; 656 else stepval = (m-1)*bs; 657 658 for (k=0; k<m; k++) { /* loop over added rows */ 659 row = im[k]; 660 if (row < 0) continue; 661 #if defined(PETSC_USE_DEBUG) 662 if (row >= a->mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,a->mbs-1); 663 #endif 664 rp = aj + ai[row]; 665 ap = aa + bs2*ai[row]; 666 rmax = imax[row]; 667 nrow = ailen[row]; 668 low = 0; 669 high = nrow; 670 for (l=0; l<n; l++) { /* loop over added columns */ 671 if (in[l] < 0) continue; 672 col = in[l]; 673 #if defined(PETSC_USE_DEBUG) 674 if (col >= a->nbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",col,a->nbs-1); 675 #endif 676 if (col < row) { 677 if (a->ignore_ltriangular) continue; /* ignore lower triangular block */ 678 else 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)"); 679 } 680 if (roworiented) value = v + k*(stepval+bs)*bs + l*bs; 681 else value = v + l*(stepval+bs)*bs + k*bs; 682 683 if (col <= lastcol) low = 0; 684 else high = nrow; 685 686 lastcol = col; 687 while (high-low > 7) { 688 t = (low+high)/2; 689 if (rp[t] > col) high = t; 690 else low = t; 691 } 692 for (i=low; i<high; i++) { 693 if (rp[i] > col) break; 694 if (rp[i] == col) { 695 bap = ap + bs2*i; 696 if (roworiented) { 697 if (is == ADD_VALUES) { 698 for (ii=0; ii<bs; ii++,value+=stepval) { 699 for (jj=ii; jj<bs2; jj+=bs) { 700 bap[jj] += *value++; 701 } 702 } 703 } else { 704 for (ii=0; ii<bs; ii++,value+=stepval) { 705 for (jj=ii; jj<bs2; jj+=bs) { 706 bap[jj] = *value++; 707 } 708 } 709 } 710 } else { 711 if (is == ADD_VALUES) { 712 for (ii=0; ii<bs; ii++,value+=stepval) { 713 for (jj=0; jj<bs; jj++) { 714 *bap++ += *value++; 715 } 716 } 717 } else { 718 for (ii=0; ii<bs; ii++,value+=stepval) { 719 for (jj=0; jj<bs; jj++) { 720 *bap++ = *value++; 721 } 722 } 723 } 724 } 725 goto noinsert2; 726 } 727 } 728 if (nonew == 1) goto noinsert2; 729 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 730 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 731 N = nrow++ - 1; high++; 732 /* shift up all the later entries in this row */ 733 for (ii=N; ii>=i; ii--) { 734 rp[ii+1] = rp[ii]; 735 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 736 } 737 if (N >= i) { 738 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 739 } 740 rp[i] = col; 741 bap = ap + bs2*i; 742 if (roworiented) { 743 for (ii=0; ii<bs; ii++,value+=stepval) { 744 for (jj=ii; jj<bs2; jj+=bs) { 745 bap[jj] = *value++; 746 } 747 } 748 } else { 749 for (ii=0; ii<bs; ii++,value+=stepval) { 750 for (jj=0; jj<bs; jj++) { 751 *bap++ = *value++; 752 } 753 } 754 } 755 noinsert2:; 756 low = i; 757 } 758 ailen[row] = nrow; 759 } 760 PetscFunctionReturn(0); 761 } 762 763 /* 764 This is not yet used 765 */ 766 #undef __FUNCT__ 767 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode" 768 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) 769 { 770 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 771 PetscErrorCode ierr; 772 const PetscInt *ai = a->i, *aj = a->j,*cols; 773 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 774 PetscBool flag; 775 776 PetscFunctionBegin; 777 ierr = PetscMalloc(m*sizeof(PetscInt),&ns);CHKERRQ(ierr); 778 while (i < m) { 779 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 780 /* Limits the number of elements in a node to 'a->inode.limit' */ 781 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 782 nzy = ai[j+1] - ai[j]; 783 if (nzy != (nzx - j + i)) break; 784 ierr = PetscMemcmp(aj + ai[i] + j - i,aj + ai[j],nzy*sizeof(PetscInt),&flag);CHKERRQ(ierr); 785 if (!flag) break; 786 } 787 ns[node_count++] = blk_size; 788 789 i = j; 790 } 791 if (!a->inode.size && m && node_count > .9*m) { 792 ierr = PetscFree(ns);CHKERRQ(ierr); 793 ierr = PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);CHKERRQ(ierr); 794 } else { 795 a->inode.node_count = node_count; 796 797 ierr = PetscMalloc(node_count*sizeof(PetscInt),&a->inode.size);CHKERRQ(ierr); 798 ierr = PetscLogObjectMemory(A,node_count*sizeof(PetscInt));CHKERRQ(ierr); 799 ierr = PetscMemcpy(a->inode.size,ns,node_count*sizeof(PetscInt));CHKERRQ(ierr); 800 ierr = PetscFree(ns);CHKERRQ(ierr); 801 ierr = PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);CHKERRQ(ierr); 802 803 /* count collections of adjacent columns in each inode */ 804 row = 0; 805 cnt = 0; 806 for (i=0; i<node_count; i++) { 807 cols = aj + ai[row] + a->inode.size[i]; 808 nz = ai[row+1] - ai[row] - a->inode.size[i]; 809 for (j=1; j<nz; j++) { 810 if (cols[j] != cols[j-1]+1) cnt++; 811 } 812 cnt++; 813 row += a->inode.size[i]; 814 } 815 ierr = PetscMalloc(2*cnt*sizeof(PetscInt),&counts);CHKERRQ(ierr); 816 cnt = 0; 817 row = 0; 818 for (i=0; i<node_count; i++) { 819 cols = aj + ai[row] + a->inode.size[i]; 820 CHKMEMQ; 821 counts[2*cnt] = cols[0]; 822 CHKMEMQ; 823 nz = ai[row+1] - ai[row] - a->inode.size[i]; 824 cnt2 = 1; 825 for (j=1; j<nz; j++) { 826 if (cols[j] != cols[j-1]+1) { 827 CHKMEMQ; 828 counts[2*(cnt++)+1] = cnt2; 829 counts[2*cnt] = cols[j]; 830 CHKMEMQ; 831 cnt2 = 1; 832 } else cnt2++; 833 } 834 CHKMEMQ; 835 counts[2*(cnt++)+1] = cnt2; 836 CHKMEMQ; 837 row += a->inode.size[i]; 838 } 839 ierr = PetscIntView(2*cnt,counts,0);CHKERRQ(ierr); 840 } 841 PetscFunctionReturn(0); 842 } 843 844 #undef __FUNCT__ 845 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 846 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 847 { 848 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 849 PetscErrorCode ierr; 850 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 851 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 852 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 853 MatScalar *aa = a->a,*ap; 854 855 PetscFunctionBegin; 856 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 857 858 if (m) rmax = ailen[0]; 859 for (i=1; i<mbs; i++) { 860 /* move each row back by the amount of empty slots (fshift) before it*/ 861 fshift += imax[i-1] - ailen[i-1]; 862 rmax = PetscMax(rmax,ailen[i]); 863 if (fshift) { 864 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 865 N = ailen[i]; 866 for (j=0; j<N; j++) { 867 ip[j-fshift] = ip[j]; 868 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 869 } 870 } 871 ai[i] = ai[i-1] + ailen[i-1]; 872 } 873 if (mbs) { 874 fshift += imax[mbs-1] - ailen[mbs-1]; 875 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 876 } 877 /* reset ilen and imax for each row */ 878 for (i=0; i<mbs; i++) { 879 ailen[i] = imax[i] = ai[i+1] - ai[i]; 880 } 881 a->nz = ai[mbs]; 882 883 /* diagonals may have moved, reset it */ 884 if (a->diag) { 885 ierr = PetscMemcpy(a->diag,ai,mbs*sizeof(PetscInt));CHKERRQ(ierr); 886 } 887 if (fshift && a->nounused == -1) 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); 888 889 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); 890 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 891 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 892 893 A->info.mallocs += a->reallocs; 894 a->reallocs = 0; 895 A->info.nz_unneeded = (PetscReal)fshift*bs2; 896 a->idiagvalid = PETSC_FALSE; 897 898 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 899 if (a->jshort && a->free_jshort) { 900 /* when matrix data structure is changed, previous jshort must be replaced */ 901 ierr = PetscFree(a->jshort);CHKERRQ(ierr); 902 } 903 ierr = PetscMalloc(a->i[A->rmap->n]*sizeof(unsigned short),&a->jshort);CHKERRQ(ierr); 904 ierr = PetscLogObjectMemory(A,a->i[A->rmap->n]*sizeof(unsigned short));CHKERRQ(ierr); 905 for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i]; 906 A->ops->mult = MatMult_SeqSBAIJ_1_ushort; 907 A->ops->sor = MatSOR_SeqSBAIJ_ushort; 908 a->free_jshort = PETSC_TRUE; 909 } 910 PetscFunctionReturn(0); 911 } 912 913 /* 914 This function returns an array of flags which indicate the locations of contiguous 915 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 916 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 917 Assume: sizes should be long enough to hold all the values. 918 */ 919 #undef __FUNCT__ 920 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 921 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 922 { 923 PetscInt i,j,k,row; 924 PetscBool flg; 925 926 PetscFunctionBegin; 927 for (i=0,j=0; i<n; j++) { 928 row = idx[i]; 929 if (row%bs!=0) { /* Not the begining of a block */ 930 sizes[j] = 1; 931 i++; 932 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 933 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 934 i++; 935 } else { /* Begining of the block, so check if the complete block exists */ 936 flg = PETSC_TRUE; 937 for (k=1; k<bs; k++) { 938 if (row+k != idx[i+k]) { /* break in the block */ 939 flg = PETSC_FALSE; 940 break; 941 } 942 } 943 if (flg) { /* No break in the bs */ 944 sizes[j] = bs; 945 i += bs; 946 } else { 947 sizes[j] = 1; 948 i++; 949 } 950 } 951 } 952 *bs_max = j; 953 PetscFunctionReturn(0); 954 } 955 956 957 /* Only add/insert a(i,j) with i<=j (blocks). 958 Any a(i,j) with i>j input by user is ingored. 959 */ 960 961 #undef __FUNCT__ 962 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 963 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 964 { 965 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 966 PetscErrorCode ierr; 967 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 968 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 969 PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 970 PetscInt ridx,cidx,bs2=a->bs2; 971 MatScalar *ap,value,*aa=a->a,*bap; 972 973 PetscFunctionBegin; 974 if (v) PetscValidScalarPointer(v,6); 975 for (k=0; k<m; k++) { /* loop over added rows */ 976 row = im[k]; /* row number */ 977 brow = row/bs; /* block row number */ 978 if (row < 0) continue; 979 #if defined(PETSC_USE_DEBUG) 980 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); 981 #endif 982 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 983 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 984 rmax = imax[brow]; /* maximum space allocated for this row */ 985 nrow = ailen[brow]; /* actual length of this row */ 986 low = 0; 987 988 for (l=0; l<n; l++) { /* loop over added columns */ 989 if (in[l] < 0) continue; 990 #if defined(PETSC_USE_DEBUG) 991 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); 992 #endif 993 col = in[l]; 994 bcol = col/bs; /* block col number */ 995 996 if (brow > bcol) { 997 if (a->ignore_ltriangular) continue; /* ignore lower triangular values */ 998 else 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)"); 999 } 1000 1001 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 1002 if ((brow==bcol && ridx<=cidx) || (brow<bcol)) { 1003 /* element value a(k,l) */ 1004 if (roworiented) value = v[l + k*n]; 1005 else value = v[k + l*m]; 1006 1007 /* move pointer bap to a(k,l) quickly and add/insert value */ 1008 if (col <= lastcol) low = 0; 1009 high = nrow; 1010 lastcol = col; 1011 while (high-low > 7) { 1012 t = (low+high)/2; 1013 if (rp[t] > bcol) high = t; 1014 else low = t; 1015 } 1016 for (i=low; i<high; i++) { 1017 if (rp[i] > bcol) break; 1018 if (rp[i] == bcol) { 1019 bap = ap + bs2*i + bs*cidx + ridx; 1020 if (is == ADD_VALUES) *bap += value; 1021 else *bap = value; 1022 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1023 if (brow == bcol && ridx < cidx) { 1024 bap = ap + bs2*i + bs*ridx + cidx; 1025 if (is == ADD_VALUES) *bap += value; 1026 else *bap = value; 1027 } 1028 goto noinsert1; 1029 } 1030 } 1031 1032 if (nonew == 1) goto noinsert1; 1033 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 1034 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 1035 1036 N = nrow++ - 1; high++; 1037 /* shift up all the later entries in this row */ 1038 for (ii=N; ii>=i; ii--) { 1039 rp[ii+1] = rp[ii]; 1040 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 1041 } 1042 if (N>=i) { 1043 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1044 } 1045 rp[i] = bcol; 1046 ap[bs2*i + bs*cidx + ridx] = value; 1047 noinsert1:; 1048 low = i; 1049 } 1050 } /* end of loop over added columns */ 1051 ailen[brow] = nrow; 1052 } /* end of loop over added rows */ 1053 PetscFunctionReturn(0); 1054 } 1055 1056 #undef __FUNCT__ 1057 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 1058 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 1059 { 1060 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1061 Mat outA; 1062 PetscErrorCode ierr; 1063 PetscBool row_identity; 1064 1065 PetscFunctionBegin; 1066 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1067 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1068 if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported"); 1069 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()! */ 1070 1071 outA = inA; 1072 inA->factortype = MAT_FACTOR_ICC; 1073 1074 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1075 ierr = MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);CHKERRQ(ierr); 1076 1077 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1078 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1079 a->row = row; 1080 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1081 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1082 a->col = row; 1083 1084 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1085 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1086 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1087 1088 if (!a->solve_work) { 1089 ierr = PetscMalloc((inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1090 ierr = PetscLogObjectMemory(inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1091 } 1092 1093 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1094 PetscFunctionReturn(0); 1095 } 1096 1097 EXTERN_C_BEGIN 1098 #undef __FUNCT__ 1099 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1100 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1101 { 1102 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data; 1103 PetscInt i,nz,n; 1104 PetscErrorCode ierr; 1105 1106 PetscFunctionBegin; 1107 nz = baij->maxnz; 1108 n = mat->cmap->n; 1109 for (i=0; i<nz; i++) baij->j[i] = indices[i]; 1110 1111 baij->nz = nz; 1112 for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i]; 1113 1114 ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1115 PetscFunctionReturn(0); 1116 } 1117 EXTERN_C_END 1118 1119 #undef __FUNCT__ 1120 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1121 /*@ 1122 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1123 in the matrix. 1124 1125 Input Parameters: 1126 + mat - the SeqSBAIJ matrix 1127 - indices - the column indices 1128 1129 Level: advanced 1130 1131 Notes: 1132 This can be called if you have precomputed the nonzero structure of the 1133 matrix and want to provide it to the matrix object to improve the performance 1134 of the MatSetValues() operation. 1135 1136 You MUST have set the correct numbers of nonzeros per row in the call to 1137 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1138 1139 MUST be called before any calls to MatSetValues() 1140 1141 .seealso: MatCreateSeqSBAIJ 1142 @*/ 1143 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1144 { 1145 PetscErrorCode ierr; 1146 1147 PetscFunctionBegin; 1148 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1149 PetscValidPointer(indices,2); 1150 ierr = PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 1151 PetscFunctionReturn(0); 1152 } 1153 1154 #undef __FUNCT__ 1155 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1156 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1157 { 1158 PetscErrorCode ierr; 1159 1160 PetscFunctionBegin; 1161 /* If the two matrices have the same copy implementation, use fast copy. */ 1162 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1163 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1164 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1165 1166 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"); 1167 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->N])*sizeof(PetscScalar));CHKERRQ(ierr); 1168 } else { 1169 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1170 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1171 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1172 } 1173 PetscFunctionReturn(0); 1174 } 1175 1176 #undef __FUNCT__ 1177 #define __FUNCT__ "MatSetUp_SeqSBAIJ" 1178 PetscErrorCode MatSetUp_SeqSBAIJ(Mat A) 1179 { 1180 PetscErrorCode ierr; 1181 1182 PetscFunctionBegin; 1183 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,A->rmap->bs,PETSC_DEFAULT,0);CHKERRQ(ierr); 1184 PetscFunctionReturn(0); 1185 } 1186 1187 #undef __FUNCT__ 1188 #define __FUNCT__ "MatSeqSBAIJGetArray_SeqSBAIJ" 1189 PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1190 { 1191 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1192 1193 PetscFunctionBegin; 1194 *array = a->a; 1195 PetscFunctionReturn(0); 1196 } 1197 1198 #undef __FUNCT__ 1199 #define __FUNCT__ "MatSeqSBAIJRestoreArray_SeqSBAIJ" 1200 PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1201 { 1202 PetscFunctionBegin; 1203 PetscFunctionReturn(0); 1204 } 1205 1206 #undef __FUNCT__ 1207 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1208 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1209 { 1210 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data; 1211 PetscErrorCode ierr; 1212 PetscInt i,bs=Y->rmap->bs,bs2=bs*bs,j; 1213 PetscBLASInt one = 1; 1214 1215 PetscFunctionBegin; 1216 if (str == SAME_NONZERO_PATTERN) { 1217 PetscScalar alpha = a; 1218 PetscBLASInt bnz; 1219 ierr = PetscBLASIntCast(x->nz*bs2,&bnz);CHKERRQ(ierr); 1220 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 1221 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1222 if (y->xtoy && y->XtoY != X) { 1223 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 1224 ierr = MatDestroy(&y->XtoY);CHKERRQ(ierr); 1225 } 1226 if (!y->xtoy) { /* get xtoy */ 1227 ierr = MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 1228 y->XtoY = X; 1229 } 1230 for (i=0; i<x->nz; i++) { 1231 j = 0; 1232 while (j < bs2) { 1233 y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]); 1234 j++; 1235 } 1236 } 1237 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); 1238 } else { 1239 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1240 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1241 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1242 } 1243 PetscFunctionReturn(0); 1244 } 1245 1246 #undef __FUNCT__ 1247 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1248 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1249 { 1250 PetscFunctionBegin; 1251 *flg = PETSC_TRUE; 1252 PetscFunctionReturn(0); 1253 } 1254 1255 #undef __FUNCT__ 1256 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1257 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg) 1258 { 1259 PetscFunctionBegin; 1260 *flg = PETSC_TRUE; 1261 PetscFunctionReturn(0); 1262 } 1263 1264 #undef __FUNCT__ 1265 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1266 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1267 { 1268 PetscFunctionBegin; 1269 *flg = PETSC_FALSE; 1270 PetscFunctionReturn(0); 1271 } 1272 1273 #undef __FUNCT__ 1274 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1275 PetscErrorCode MatRealPart_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] = PetscRealPart(aa[i]); 1283 PetscFunctionReturn(0); 1284 } 1285 1286 #undef __FUNCT__ 1287 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1288 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1289 { 1290 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1291 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1292 MatScalar *aa = a->a; 1293 1294 PetscFunctionBegin; 1295 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1296 PetscFunctionReturn(0); 1297 } 1298 1299 #undef __FUNCT__ 1300 #define __FUNCT__ "MatZeroRowsColumns_SeqSBAIJ" 1301 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b) 1302 { 1303 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data; 1304 PetscErrorCode ierr; 1305 PetscInt i,j,k,count; 1306 PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col; 1307 PetscScalar zero = 0.0; 1308 MatScalar *aa; 1309 const PetscScalar *xx; 1310 PetscScalar *bb; 1311 PetscBool *zeroed,vecs = PETSC_FALSE; 1312 1313 PetscFunctionBegin; 1314 /* fix right hand side if needed */ 1315 if (x && b) { 1316 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1317 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1318 vecs = PETSC_TRUE; 1319 } 1320 A->same_nonzero = PETSC_TRUE; 1321 1322 /* zero the columns */ 1323 ierr = PetscMalloc(A->rmap->n*sizeof(PetscBool),&zeroed);CHKERRQ(ierr); 1324 ierr = PetscMemzero(zeroed,A->rmap->n*sizeof(PetscBool));CHKERRQ(ierr); 1325 for (i=0; i<is_n; i++) { 1326 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]); 1327 zeroed[is_idx[i]] = PETSC_TRUE; 1328 } 1329 if (vecs) { 1330 for (i=0; i<A->rmap->N; i++) { 1331 row = i/bs; 1332 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1333 for (k=0; k<bs; k++) { 1334 col = bs*baij->j[j] + k; 1335 if (col <= i) continue; 1336 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1337 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col]; 1338 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i]; 1339 } 1340 } 1341 } 1342 for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]]; 1343 } 1344 1345 for (i=0; i<A->rmap->N; i++) { 1346 if (!zeroed[i]) { 1347 row = i/bs; 1348 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1349 for (k=0; k<bs; k++) { 1350 col = bs*baij->j[j] + k; 1351 if (zeroed[col]) { 1352 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1353 aa[0] = 0.0; 1354 } 1355 } 1356 } 1357 } 1358 } 1359 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1360 if (vecs) { 1361 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1362 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1363 } 1364 1365 /* zero the rows */ 1366 for (i=0; i<is_n; i++) { 1367 row = is_idx[i]; 1368 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1369 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1370 for (k=0; k<count; k++) { 1371 aa[0] = zero; 1372 aa += bs; 1373 } 1374 if (diag != 0.0) { 1375 ierr = (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr); 1376 } 1377 } 1378 ierr = MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1379 PetscFunctionReturn(0); 1380 } 1381 1382 /* -------------------------------------------------------------------*/ 1383 static struct _MatOps MatOps_Values = { 1384 MatSetValues_SeqSBAIJ, 1385 MatGetRow_SeqSBAIJ, 1386 MatRestoreRow_SeqSBAIJ, 1387 MatMult_SeqSBAIJ_N, 1388 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1389 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1390 MatMultAdd_SeqSBAIJ_N, 1391 0, 1392 0, 1393 0, 1394 /* 10*/ 0, 1395 0, 1396 MatCholeskyFactor_SeqSBAIJ, 1397 MatSOR_SeqSBAIJ, 1398 MatTranspose_SeqSBAIJ, 1399 /* 15*/ MatGetInfo_SeqSBAIJ, 1400 MatEqual_SeqSBAIJ, 1401 MatGetDiagonal_SeqSBAIJ, 1402 MatDiagonalScale_SeqSBAIJ, 1403 MatNorm_SeqSBAIJ, 1404 /* 20*/ 0, 1405 MatAssemblyEnd_SeqSBAIJ, 1406 MatSetOption_SeqSBAIJ, 1407 MatZeroEntries_SeqSBAIJ, 1408 /* 24*/ 0, 1409 0, 1410 0, 1411 0, 1412 0, 1413 /* 29*/ MatSetUp_SeqSBAIJ, 1414 0, 1415 0, 1416 0, 1417 0, 1418 /* 34*/ MatDuplicate_SeqSBAIJ, 1419 0, 1420 0, 1421 0, 1422 MatICCFactor_SeqSBAIJ, 1423 /* 39*/ MatAXPY_SeqSBAIJ, 1424 MatGetSubMatrices_SeqSBAIJ, 1425 MatIncreaseOverlap_SeqSBAIJ, 1426 MatGetValues_SeqSBAIJ, 1427 MatCopy_SeqSBAIJ, 1428 /* 44*/ 0, 1429 MatScale_SeqSBAIJ, 1430 0, 1431 0, 1432 MatZeroRowsColumns_SeqSBAIJ, 1433 /* 49*/ 0, 1434 MatGetRowIJ_SeqSBAIJ, 1435 MatRestoreRowIJ_SeqSBAIJ, 1436 0, 1437 0, 1438 /* 54*/ 0, 1439 0, 1440 0, 1441 0, 1442 MatSetValuesBlocked_SeqSBAIJ, 1443 /* 59*/ MatGetSubMatrix_SeqSBAIJ, 1444 0, 1445 0, 1446 0, 1447 0, 1448 /* 64*/ 0, 1449 0, 1450 0, 1451 0, 1452 0, 1453 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1454 0, 1455 0, 1456 0, 1457 0, 1458 /* 74*/ 0, 1459 0, 1460 0, 1461 0, 1462 0, 1463 /* 79*/ 0, 1464 0, 1465 0, 1466 MatGetInertia_SeqSBAIJ, 1467 MatLoad_SeqSBAIJ, 1468 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1469 MatIsHermitian_SeqSBAIJ, 1470 MatIsStructurallySymmetric_SeqSBAIJ, 1471 0, 1472 0, 1473 /* 89*/ 0, 1474 0, 1475 0, 1476 0, 1477 0, 1478 /* 94*/ 0, 1479 0, 1480 0, 1481 0, 1482 0, 1483 /* 99*/ 0, 1484 0, 1485 0, 1486 0, 1487 0, 1488 /*104*/ 0, 1489 MatRealPart_SeqSBAIJ, 1490 MatImaginaryPart_SeqSBAIJ, 1491 MatGetRowUpperTriangular_SeqSBAIJ, 1492 MatRestoreRowUpperTriangular_SeqSBAIJ, 1493 /*109*/ 0, 1494 0, 1495 0, 1496 0, 1497 MatMissingDiagonal_SeqSBAIJ, 1498 /*114*/ 0, 1499 0, 1500 0, 1501 0, 1502 0, 1503 /*119*/ 0, 1504 0, 1505 0, 1506 0 1507 }; 1508 1509 EXTERN_C_BEGIN 1510 #undef __FUNCT__ 1511 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1512 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1513 { 1514 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1515 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1516 PetscErrorCode ierr; 1517 1518 PetscFunctionBegin; 1519 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1520 1521 /* allocate space for values if not already there */ 1522 if (!aij->saved_values) { 1523 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1524 } 1525 1526 /* copy values over */ 1527 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1528 PetscFunctionReturn(0); 1529 } 1530 EXTERN_C_END 1531 1532 EXTERN_C_BEGIN 1533 #undef __FUNCT__ 1534 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1535 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1536 { 1537 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1538 PetscErrorCode ierr; 1539 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1540 1541 PetscFunctionBegin; 1542 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1543 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1544 1545 /* copy values over */ 1546 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1547 PetscFunctionReturn(0); 1548 } 1549 EXTERN_C_END 1550 1551 EXTERN_C_BEGIN 1552 #undef __FUNCT__ 1553 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1554 PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1555 { 1556 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1557 PetscErrorCode ierr; 1558 PetscInt i,mbs,bs2; 1559 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1560 1561 PetscFunctionBegin; 1562 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1563 B->preallocated = PETSC_TRUE; 1564 1565 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1566 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1567 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1568 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1569 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1570 1571 mbs = B->rmap->N/bs; 1572 bs2 = bs*bs; 1573 1574 if (mbs*bs != B->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1575 1576 if (nz == MAT_SKIP_ALLOCATION) { 1577 skipallocation = PETSC_TRUE; 1578 nz = 0; 1579 } 1580 1581 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1582 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1583 if (nnz) { 1584 for (i=0; i<mbs; i++) { 1585 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]); 1586 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); 1587 } 1588 } 1589 1590 B->ops->mult = MatMult_SeqSBAIJ_N; 1591 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1592 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1593 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1594 1595 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1596 if (!flg) { 1597 switch (bs) { 1598 case 1: 1599 B->ops->mult = MatMult_SeqSBAIJ_1; 1600 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1601 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1602 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1603 break; 1604 case 2: 1605 B->ops->mult = MatMult_SeqSBAIJ_2; 1606 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1607 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1608 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1609 break; 1610 case 3: 1611 B->ops->mult = MatMult_SeqSBAIJ_3; 1612 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1613 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1614 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1615 break; 1616 case 4: 1617 B->ops->mult = MatMult_SeqSBAIJ_4; 1618 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1619 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1620 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1621 break; 1622 case 5: 1623 B->ops->mult = MatMult_SeqSBAIJ_5; 1624 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1625 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1626 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1627 break; 1628 case 6: 1629 B->ops->mult = MatMult_SeqSBAIJ_6; 1630 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1631 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1632 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1633 break; 1634 case 7: 1635 B->ops->mult = MatMult_SeqSBAIJ_7; 1636 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1637 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1638 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1639 break; 1640 } 1641 } 1642 1643 b->mbs = mbs; 1644 b->nbs = mbs; 1645 if (!skipallocation) { 1646 if (!b->imax) { 1647 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1648 1649 b->free_imax_ilen = PETSC_TRUE; 1650 1651 ierr = PetscLogObjectMemory(B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1652 } 1653 if (!nnz) { 1654 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1655 else if (nz <= 0) nz = 1; 1656 for (i=0; i<mbs; i++) b->imax[i] = nz; 1657 nz = nz*mbs; /* total nz */ 1658 } else { 1659 nz = 0; 1660 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1661 } 1662 /* b->ilen will count nonzeros in each block row so far. */ 1663 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1664 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1665 1666 /* allocate the matrix space */ 1667 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1668 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->N+1,PetscInt,&b->i);CHKERRQ(ierr); 1669 ierr = PetscLogObjectMemory(B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1670 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1671 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1672 1673 b->singlemalloc = PETSC_TRUE; 1674 1675 /* pointer to beginning of each row */ 1676 b->i[0] = 0; 1677 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1678 1679 b->free_a = PETSC_TRUE; 1680 b->free_ij = PETSC_TRUE; 1681 } else { 1682 b->free_a = PETSC_FALSE; 1683 b->free_ij = PETSC_FALSE; 1684 } 1685 1686 B->rmap->bs = bs; 1687 b->bs2 = bs2; 1688 b->nz = 0; 1689 b->maxnz = nz; 1690 1691 b->inew = 0; 1692 b->jnew = 0; 1693 b->anew = 0; 1694 b->a2anew = 0; 1695 b->permute = PETSC_FALSE; 1696 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 1697 PetscFunctionReturn(0); 1698 } 1699 EXTERN_C_END 1700 1701 /* 1702 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1703 */ 1704 #undef __FUNCT__ 1705 #define __FUNCT__ "MatSeqSBAIJSetNumericFactorization_inplace" 1706 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1707 { 1708 PetscErrorCode ierr; 1709 PetscBool flg = PETSC_FALSE; 1710 PetscInt bs = B->rmap->bs; 1711 1712 PetscFunctionBegin; 1713 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1714 if (flg) bs = 8; 1715 1716 if (!natural) { 1717 switch (bs) { 1718 case 1: 1719 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1720 break; 1721 case 2: 1722 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1723 break; 1724 case 3: 1725 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1726 break; 1727 case 4: 1728 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1729 break; 1730 case 5: 1731 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1732 break; 1733 case 6: 1734 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1735 break; 1736 case 7: 1737 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1738 break; 1739 default: 1740 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1741 break; 1742 } 1743 } else { 1744 switch (bs) { 1745 case 1: 1746 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1747 break; 1748 case 2: 1749 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1750 break; 1751 case 3: 1752 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1753 break; 1754 case 4: 1755 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1756 break; 1757 case 5: 1758 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1759 break; 1760 case 6: 1761 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1762 break; 1763 case 7: 1764 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1765 break; 1766 default: 1767 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1768 break; 1769 } 1770 } 1771 PetscFunctionReturn(0); 1772 } 1773 1774 EXTERN_C_BEGIN 1775 extern PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1776 extern PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1777 EXTERN_C_END 1778 1779 1780 EXTERN_C_BEGIN 1781 #undef __FUNCT__ 1782 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1783 PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1784 { 1785 PetscInt n = A->rmap->n; 1786 PetscErrorCode ierr; 1787 1788 PetscFunctionBegin; 1789 #if defined(PETSC_USE_COMPLEX) 1790 if (A->hermitian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported"); 1791 #endif 1792 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1793 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1794 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1795 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1796 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1797 1798 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1799 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1800 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1801 (*B)->factortype = ftype; 1802 PetscFunctionReturn(0); 1803 } 1804 EXTERN_C_END 1805 1806 EXTERN_C_BEGIN 1807 #undef __FUNCT__ 1808 #define __FUNCT__ "MatGetFactorAvailable_seqsbaij_petsc" 1809 PetscErrorCode MatGetFactorAvailable_seqsbaij_petsc(Mat A,MatFactorType ftype,PetscBool *flg) 1810 { 1811 PetscFunctionBegin; 1812 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1813 *flg = PETSC_TRUE; 1814 } else { 1815 *flg = PETSC_FALSE; 1816 } 1817 PetscFunctionReturn(0); 1818 } 1819 EXTERN_C_END 1820 1821 EXTERN_C_BEGIN 1822 #if defined(PETSC_HAVE_MUMPS) 1823 extern PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*); 1824 #endif 1825 #if defined(PETSC_HAVE_PASTIX) 1826 extern PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat,MatFactorType,Mat*); 1827 #endif 1828 #if defined(PETSC_HAVE_CHOLMOD) 1829 extern PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat,MatFactorType,Mat*); 1830 #endif 1831 extern PetscErrorCode MatGetFactor_seqsbaij_sbstrm(Mat,MatFactorType,Mat*); 1832 EXTERN_C_END 1833 1834 /*MC 1835 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1836 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1837 1838 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1839 can call MatSetOption(Mat, MAT_HERMITIAN); after MatAssemblyEnd() 1840 1841 Options Database Keys: 1842 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1843 1844 Notes: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1845 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1846 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. 1847 1848 1849 Level: beginner 1850 1851 .seealso: MatCreateSeqSBAIJ 1852 M*/ 1853 1854 EXTERN_C_BEGIN 1855 extern PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat, MatType,MatReuse,Mat*); 1856 EXTERN_C_END 1857 1858 1859 EXTERN_C_BEGIN 1860 #undef __FUNCT__ 1861 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1862 PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1863 { 1864 Mat_SeqSBAIJ *b; 1865 PetscErrorCode ierr; 1866 PetscMPIInt size; 1867 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1868 1869 PetscFunctionBegin; 1870 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1871 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1872 1873 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1874 B->data = (void*)b; 1875 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1876 1877 B->ops->destroy = MatDestroy_SeqSBAIJ; 1878 B->ops->view = MatView_SeqSBAIJ; 1879 b->row = 0; 1880 b->icol = 0; 1881 b->reallocs = 0; 1882 b->saved_values = 0; 1883 b->inode.limit = 5; 1884 b->inode.max_limit = 5; 1885 1886 b->roworiented = PETSC_TRUE; 1887 b->nonew = 0; 1888 b->diag = 0; 1889 b->solve_work = 0; 1890 b->mult_work = 0; 1891 B->spptr = 0; 1892 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1893 b->keepnonzeropattern = PETSC_FALSE; 1894 b->xtoy = 0; 1895 b->XtoY = 0; 1896 1897 b->inew = 0; 1898 b->jnew = 0; 1899 b->anew = 0; 1900 b->a2anew = 0; 1901 b->permute = PETSC_FALSE; 1902 1903 b->ignore_ltriangular = PETSC_TRUE; 1904 1905 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,PETSC_NULL);CHKERRQ(ierr); 1906 1907 b->getrow_utriangular = PETSC_FALSE; 1908 1909 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,PETSC_NULL);CHKERRQ(ierr); 1910 1911 #if defined(PETSC_HAVE_PASTIX) 1912 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 1913 "MatGetFactor_seqsbaij_pastix", 1914 MatGetFactor_seqsbaij_pastix);CHKERRQ(ierr); 1915 #endif 1916 #if defined(PETSC_HAVE_MUMPS) 1917 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 1918 "MatGetFactor_sbaij_mumps", 1919 MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 1920 #endif 1921 #if defined(PETSC_HAVE_CHOLMOD) 1922 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 1923 "MatGetFactor_seqsbaij_cholmod", 1924 MatGetFactor_seqsbaij_cholmod);CHKERRQ(ierr); 1925 #endif 1926 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 1927 "MatGetFactorAvailable_seqsbaij_petsc", 1928 MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1929 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 1930 "MatGetFactor_seqsbaij_petsc", 1931 MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1932 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_sbstrm_C", 1933 "MatGetFactor_seqsbaij_sbstrm", 1934 MatGetFactor_seqsbaij_sbstrm);CHKERRQ(ierr); 1935 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1936 "MatStoreValues_SeqSBAIJ", 1937 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1938 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1939 "MatRetrieveValues_SeqSBAIJ", 1940 MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1941 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1942 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1943 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1944 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1945 "MatConvert_SeqSBAIJ_SeqAIJ", 1946 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1947 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1948 "MatConvert_SeqSBAIJ_SeqBAIJ", 1949 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1950 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1951 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1952 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1953 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqsbstrm_C", 1954 "MatConvert_SeqSBAIJ_SeqSBSTRM", 1955 MatConvert_SeqSBAIJ_SeqSBSTRM);CHKERRQ(ierr); 1956 1957 B->symmetric = PETSC_TRUE; 1958 B->structurally_symmetric = PETSC_TRUE; 1959 B->symmetric_set = PETSC_TRUE; 1960 B->structurally_symmetric_set = PETSC_TRUE; 1961 1962 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1963 1964 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1965 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",PETSC_NULL,no_unroll,&no_unroll,PETSC_NULL);CHKERRQ(ierr); 1966 if (no_unroll) { 1967 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr); 1968 } 1969 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",PETSC_NULL,no_inode,&no_inode,PETSC_NULL);CHKERRQ(ierr); 1970 if (no_inode) { 1971 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr); 1972 } 1973 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); 1974 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1975 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 1976 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1977 PetscFunctionReturn(0); 1978 } 1979 EXTERN_C_END 1980 1981 #undef __FUNCT__ 1982 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1983 /*@C 1984 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1985 compressed row) format. For good matrix assembly performance the 1986 user should preallocate the matrix storage by setting the parameter nz 1987 (or the array nnz). By setting these parameters accurately, performance 1988 during matrix assembly can be increased by more than a factor of 50. 1989 1990 Collective on Mat 1991 1992 Input Parameters: 1993 + A - the symmetric matrix 1994 . bs - size of block 1995 . nz - number of block nonzeros per block row (same for all rows) 1996 - nnz - array containing the number of block nonzeros in the upper triangular plus 1997 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1998 1999 Options Database Keys: 2000 . -mat_no_unroll - uses code that does not unroll the loops in the 2001 block calculations (much slower) 2002 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2003 2004 Level: intermediate 2005 2006 Notes: 2007 Specify the preallocated storage with either nz or nnz (not both). 2008 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2009 allocation. See the <a href="../../docs/manual.pdf#nameddest=ch_mat">Mat chapter of the users manual</a> for details. 2010 2011 You can call MatGetInfo() to get information on how effective the preallocation was; 2012 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2013 You can also run with the option -info and look for messages with the string 2014 malloc in them to see if additional memory allocation was needed. 2015 2016 If the nnz parameter is given then the nz parameter is ignored 2017 2018 2019 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2020 @*/ 2021 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2022 { 2023 PetscErrorCode ierr; 2024 2025 PetscFunctionBegin; 2026 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2027 PetscValidType(B,1); 2028 PetscValidLogicalCollectiveInt(B,bs,2); 2029 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2030 PetscFunctionReturn(0); 2031 } 2032 2033 #undef __FUNCT__ 2034 #define __FUNCT__ "MatCreateSeqSBAIJ" 2035 /*@C 2036 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2037 compressed row) format. For good matrix assembly performance the 2038 user should preallocate the matrix storage by setting the parameter nz 2039 (or the array nnz). By setting these parameters accurately, performance 2040 during matrix assembly can be increased by more than a factor of 50. 2041 2042 Collective on MPI_Comm 2043 2044 Input Parameters: 2045 + comm - MPI communicator, set to PETSC_COMM_SELF 2046 . bs - size of block 2047 . m - number of rows, or number of columns 2048 . nz - number of block nonzeros per block row (same for all rows) 2049 - nnz - array containing the number of block nonzeros in the upper triangular plus 2050 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 2051 2052 Output Parameter: 2053 . A - the symmetric matrix 2054 2055 Options Database Keys: 2056 . -mat_no_unroll - uses code that does not unroll the loops in the 2057 block calculations (much slower) 2058 . -mat_block_size - size of the blocks to use 2059 2060 Level: intermediate 2061 2062 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2063 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2064 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2065 2066 Notes: 2067 The number of rows and columns must be divisible by blocksize. 2068 This matrix type does not support complex Hermitian operation. 2069 2070 Specify the preallocated storage with either nz or nnz (not both). 2071 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2072 allocation. See the <a href="../../docs/manual.pdf#nameddest=ch_mat">Mat chapter of the users manual</a> for details. 2073 2074 If the nnz parameter is given then the nz parameter is ignored 2075 2076 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2077 @*/ 2078 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2079 { 2080 PetscErrorCode ierr; 2081 2082 PetscFunctionBegin; 2083 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2084 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2085 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2086 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2087 PetscFunctionReturn(0); 2088 } 2089 2090 #undef __FUNCT__ 2091 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 2092 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2093 { 2094 Mat C; 2095 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2096 PetscErrorCode ierr; 2097 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2098 2099 PetscFunctionBegin; 2100 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2101 2102 *B = 0; 2103 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 2104 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2105 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2106 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2107 c = (Mat_SeqSBAIJ*)C->data; 2108 2109 C->preallocated = PETSC_TRUE; 2110 C->factortype = A->factortype; 2111 c->row = 0; 2112 c->icol = 0; 2113 c->saved_values = 0; 2114 c->keepnonzeropattern = a->keepnonzeropattern; 2115 C->assembled = PETSC_TRUE; 2116 2117 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2118 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2119 c->bs2 = a->bs2; 2120 c->mbs = a->mbs; 2121 c->nbs = a->nbs; 2122 2123 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2124 c->imax = a->imax; 2125 c->ilen = a->ilen; 2126 c->free_imax_ilen = PETSC_FALSE; 2127 } else { 2128 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 2129 ierr = PetscLogObjectMemory(C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2130 for (i=0; i<mbs; i++) { 2131 c->imax[i] = a->imax[i]; 2132 c->ilen[i] = a->ilen[i]; 2133 } 2134 c->free_imax_ilen = PETSC_TRUE; 2135 } 2136 2137 /* allocate the matrix space */ 2138 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2139 ierr = PetscMalloc(bs2*nz*sizeof(MatScalar),&c->a);CHKERRQ(ierr); 2140 ierr = PetscLogObjectMemory(C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2141 c->i = a->i; 2142 c->j = a->j; 2143 c->singlemalloc = PETSC_FALSE; 2144 c->free_a = PETSC_TRUE; 2145 c->free_ij = PETSC_FALSE; 2146 c->parent = A; 2147 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2148 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2149 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2150 } else { 2151 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 2152 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2153 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2154 c->singlemalloc = PETSC_TRUE; 2155 c->free_a = PETSC_TRUE; 2156 c->free_ij = PETSC_TRUE; 2157 } 2158 if (mbs > 0) { 2159 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2160 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2161 } 2162 if (cpvalues == MAT_COPY_VALUES) { 2163 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2164 } else { 2165 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2166 } 2167 if (a->jshort) { 2168 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2169 /* if the parent matrix is reassembled, this child matrix will never notice */ 2170 ierr = PetscMalloc(nz*sizeof(unsigned short),&c->jshort);CHKERRQ(ierr); 2171 ierr = PetscLogObjectMemory(C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2172 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2173 2174 c->free_jshort = PETSC_TRUE; 2175 } 2176 } 2177 2178 c->roworiented = a->roworiented; 2179 c->nonew = a->nonew; 2180 2181 if (a->diag) { 2182 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2183 c->diag = a->diag; 2184 c->free_diag = PETSC_FALSE; 2185 } else { 2186 ierr = PetscMalloc(mbs*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2187 ierr = PetscLogObjectMemory(C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2188 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2189 c->free_diag = PETSC_TRUE; 2190 } 2191 } 2192 c->nz = a->nz; 2193 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2194 c->solve_work = 0; 2195 c->mult_work = 0; 2196 2197 *B = C; 2198 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2199 PetscFunctionReturn(0); 2200 } 2201 2202 #undef __FUNCT__ 2203 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2204 PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer) 2205 { 2206 Mat_SeqSBAIJ *a; 2207 PetscErrorCode ierr; 2208 int fd; 2209 PetscMPIInt size; 2210 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2211 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2212 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2213 PetscInt *masked,nmask,tmp,bs2,ishift; 2214 PetscScalar *aa; 2215 MPI_Comm comm = ((PetscObject)viewer)->comm; 2216 2217 PetscFunctionBegin; 2218 ierr = PetscOptionsGetInt(((PetscObject)newmat)->prefix,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 2219 bs2 = bs*bs; 2220 2221 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2222 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2223 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2224 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2225 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2226 M = header[1]; N = header[2]; nz = header[3]; 2227 2228 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2229 2230 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2231 2232 /* 2233 This code adds extra rows to make sure the number of rows is 2234 divisible by the blocksize 2235 */ 2236 mbs = M/bs; 2237 extra_rows = bs - M + bs*(mbs); 2238 if (extra_rows == bs) extra_rows = 0; 2239 else mbs++; 2240 if (extra_rows) { 2241 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2242 } 2243 2244 /* Set global sizes if not already set */ 2245 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2246 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2247 } else { /* Check if the matrix global sizes are correct */ 2248 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2249 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); 2250 } 2251 2252 /* read in row lengths */ 2253 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2254 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2255 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2256 2257 /* read in column indices */ 2258 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2259 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2260 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2261 2262 /* loop over row lengths determining block row lengths */ 2263 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2264 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2265 ierr = PetscMalloc2(mbs,PetscInt,&mask,mbs,PetscInt,&masked);CHKERRQ(ierr); 2266 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2267 rowcount = 0; 2268 nzcount = 0; 2269 for (i=0; i<mbs; i++) { 2270 nmask = 0; 2271 for (j=0; j<bs; j++) { 2272 kmax = rowlengths[rowcount]; 2273 for (k=0; k<kmax; k++) { 2274 tmp = jj[nzcount++]/bs; /* block col. index */ 2275 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2276 } 2277 rowcount++; 2278 } 2279 s_browlengths[i] += nmask; 2280 2281 /* zero out the mask elements we set */ 2282 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2283 } 2284 2285 /* Do preallocation */ 2286 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2287 a = (Mat_SeqSBAIJ*)newmat->data; 2288 2289 /* set matrix "i" values */ 2290 a->i[0] = 0; 2291 for (i=1; i<= mbs; i++) { 2292 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2293 a->ilen[i-1] = s_browlengths[i-1]; 2294 } 2295 a->nz = a->i[mbs]; 2296 2297 /* read in nonzero values */ 2298 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2299 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2300 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2301 2302 /* set "a" and "j" values into matrix */ 2303 nzcount = 0; jcount = 0; 2304 for (i=0; i<mbs; i++) { 2305 nzcountb = nzcount; 2306 nmask = 0; 2307 for (j=0; j<bs; j++) { 2308 kmax = rowlengths[i*bs+j]; 2309 for (k=0; k<kmax; k++) { 2310 tmp = jj[nzcount++]/bs; /* block col. index */ 2311 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2312 } 2313 } 2314 /* sort the masked values */ 2315 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2316 2317 /* set "j" values into matrix */ 2318 maskcount = 1; 2319 for (j=0; j<nmask; j++) { 2320 a->j[jcount++] = masked[j]; 2321 mask[masked[j]] = maskcount++; 2322 } 2323 2324 /* set "a" values into matrix */ 2325 ishift = bs2*a->i[i]; 2326 for (j=0; j<bs; j++) { 2327 kmax = rowlengths[i*bs+j]; 2328 for (k=0; k<kmax; k++) { 2329 tmp = jj[nzcountb]/bs; /* block col. index */ 2330 if (tmp >= i) { 2331 block = mask[tmp] - 1; 2332 point = jj[nzcountb] - bs*tmp; 2333 idx = ishift + bs2*block + j + bs*point; 2334 a->a[idx] = aa[nzcountb]; 2335 } 2336 nzcountb++; 2337 } 2338 } 2339 /* zero out the mask elements we set */ 2340 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2341 } 2342 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2343 2344 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2345 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2346 ierr = PetscFree(aa);CHKERRQ(ierr); 2347 ierr = PetscFree(jj);CHKERRQ(ierr); 2348 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2349 2350 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2351 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2352 PetscFunctionReturn(0); 2353 } 2354 2355 #undef __FUNCT__ 2356 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2357 /*@ 2358 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2359 (upper triangular entries in CSR format) provided by the user. 2360 2361 Collective on MPI_Comm 2362 2363 Input Parameters: 2364 + comm - must be an MPI communicator of size 1 2365 . bs - size of block 2366 . m - number of rows 2367 . n - number of columns 2368 . i - row indices 2369 . j - column indices 2370 - a - matrix values 2371 2372 Output Parameter: 2373 . mat - the matrix 2374 2375 Level: advanced 2376 2377 Notes: 2378 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2379 once the matrix is destroyed 2380 2381 You cannot set new nonzero locations into this matrix, that will generate an error. 2382 2383 The i and j indices are 0 based 2384 2385 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 2386 it is the regular CSR format excluding the lower triangular elements. 2387 2388 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ() 2389 2390 @*/ 2391 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 2392 { 2393 PetscErrorCode ierr; 2394 PetscInt ii; 2395 Mat_SeqSBAIJ *sbaij; 2396 2397 PetscFunctionBegin; 2398 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2399 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2400 2401 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2402 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2403 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2404 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2405 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2406 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2407 ierr = PetscLogObjectMemory(*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2408 2409 sbaij->i = i; 2410 sbaij->j = j; 2411 sbaij->a = a; 2412 2413 sbaij->singlemalloc = PETSC_FALSE; 2414 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2415 sbaij->free_a = PETSC_FALSE; 2416 sbaij->free_ij = PETSC_FALSE; 2417 2418 for (ii=0; ii<m; ii++) { 2419 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2420 #if defined(PETSC_USE_DEBUG) 2421 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]); 2422 #endif 2423 } 2424 #if defined(PETSC_USE_DEBUG) 2425 for (ii=0; ii<sbaij->i[m]; ii++) { 2426 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2427 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]); 2428 } 2429 #endif 2430 2431 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2432 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2433 PetscFunctionReturn(0); 2434 } 2435 2436 2437 2438 2439 2440