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