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