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