1 /* 2 Defines the basic matrix operations for the AIJ (compressed row) 3 matrix storage format. 4 */ 5 6 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 #include <petscbt.h> 9 #include <petsc/private/kernels/blocktranspose.h> 10 11 PetscErrorCode MatSeqAIJSetTypeFromOptions(Mat A) 12 { 13 PetscBool flg; 14 char type[256]; 15 16 PetscFunctionBegin; 17 PetscObjectOptionsBegin((PetscObject)A); 18 PetscCall(PetscOptionsFList("-mat_seqaij_type","Matrix SeqAIJ type","MatSeqAIJSetType",MatSeqAIJList,"seqaij",type,256,&flg)); 19 if (flg) PetscCall(MatSeqAIJSetType(A,type)); 20 PetscOptionsEnd(); 21 PetscFunctionReturn(0); 22 } 23 24 PetscErrorCode MatGetColumnReductions_SeqAIJ(Mat A,PetscInt type,PetscReal *reductions) 25 { 26 PetscInt i,m,n; 27 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 28 29 PetscFunctionBegin; 30 PetscCall(MatGetSize(A,&m,&n)); 31 PetscCall(PetscArrayzero(reductions,n)); 32 if (type == NORM_2) { 33 for (i=0; i<aij->i[m]; i++) { 34 reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 35 } 36 } else if (type == NORM_1) { 37 for (i=0; i<aij->i[m]; i++) { 38 reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]); 39 } 40 } else if (type == NORM_INFINITY) { 41 for (i=0; i<aij->i[m]; i++) { 42 reductions[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),reductions[aij->j[i]]); 43 } 44 } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) { 45 for (i=0; i<aij->i[m]; i++) { 46 reductions[aij->j[i]] += PetscRealPart(aij->a[i]); 47 } 48 } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) { 49 for (i=0; i<aij->i[m]; i++) { 50 reductions[aij->j[i]] += PetscImaginaryPart(aij->a[i]); 51 } 52 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown reduction type"); 53 54 if (type == NORM_2) { 55 for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]); 56 } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) { 57 for (i=0; i<n; i++) reductions[i] /= m; 58 } 59 PetscFunctionReturn(0); 60 } 61 62 PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A,IS *is) 63 { 64 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 65 PetscInt i,m=A->rmap->n,cnt = 0, bs = A->rmap->bs; 66 const PetscInt *jj = a->j,*ii = a->i; 67 PetscInt *rows; 68 69 PetscFunctionBegin; 70 for (i=0; i<m; i++) { 71 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 72 cnt++; 73 } 74 } 75 PetscCall(PetscMalloc1(cnt,&rows)); 76 cnt = 0; 77 for (i=0; i<m; i++) { 78 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 79 rows[cnt] = i; 80 cnt++; 81 } 82 } 83 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,is)); 84 PetscFunctionReturn(0); 85 } 86 87 PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows) 88 { 89 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 90 const MatScalar *aa; 91 PetscInt i,m=A->rmap->n,cnt = 0; 92 const PetscInt *ii = a->i,*jj = a->j,*diag; 93 PetscInt *rows; 94 95 PetscFunctionBegin; 96 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 97 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 98 diag = a->diag; 99 for (i=0; i<m; i++) { 100 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 101 cnt++; 102 } 103 } 104 PetscCall(PetscMalloc1(cnt,&rows)); 105 cnt = 0; 106 for (i=0; i<m; i++) { 107 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 108 rows[cnt++] = i; 109 } 110 } 111 *nrows = cnt; 112 *zrows = rows; 113 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 114 PetscFunctionReturn(0); 115 } 116 117 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows) 118 { 119 PetscInt nrows,*rows; 120 121 PetscFunctionBegin; 122 *zrows = NULL; 123 PetscCall(MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows)); 124 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows)); 125 PetscFunctionReturn(0); 126 } 127 128 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows) 129 { 130 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 131 const MatScalar *aa; 132 PetscInt m=A->rmap->n,cnt = 0; 133 const PetscInt *ii; 134 PetscInt n,i,j,*rows; 135 136 PetscFunctionBegin; 137 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 138 *keptrows = NULL; 139 ii = a->i; 140 for (i=0; i<m; i++) { 141 n = ii[i+1] - ii[i]; 142 if (!n) { 143 cnt++; 144 goto ok1; 145 } 146 for (j=ii[i]; j<ii[i+1]; j++) { 147 if (aa[j] != 0.0) goto ok1; 148 } 149 cnt++; 150 ok1:; 151 } 152 if (!cnt) { 153 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 154 PetscFunctionReturn(0); 155 } 156 PetscCall(PetscMalloc1(A->rmap->n-cnt,&rows)); 157 cnt = 0; 158 for (i=0; i<m; i++) { 159 n = ii[i+1] - ii[i]; 160 if (!n) continue; 161 for (j=ii[i]; j<ii[i+1]; j++) { 162 if (aa[j] != 0.0) { 163 rows[cnt++] = i; 164 break; 165 } 166 } 167 } 168 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 169 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows)); 170 PetscFunctionReturn(0); 171 } 172 173 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 174 { 175 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 176 PetscInt i,m = Y->rmap->n; 177 const PetscInt *diag; 178 MatScalar *aa; 179 const PetscScalar *v; 180 PetscBool missing; 181 182 PetscFunctionBegin; 183 if (Y->assembled) { 184 PetscCall(MatMissingDiagonal_SeqAIJ(Y,&missing,NULL)); 185 if (!missing) { 186 diag = aij->diag; 187 PetscCall(VecGetArrayRead(D,&v)); 188 PetscCall(MatSeqAIJGetArray(Y,&aa)); 189 if (is == INSERT_VALUES) { 190 for (i=0; i<m; i++) { 191 aa[diag[i]] = v[i]; 192 } 193 } else { 194 for (i=0; i<m; i++) { 195 aa[diag[i]] += v[i]; 196 } 197 } 198 PetscCall(MatSeqAIJRestoreArray(Y,&aa)); 199 PetscCall(VecRestoreArrayRead(D,&v)); 200 PetscFunctionReturn(0); 201 } 202 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 203 } 204 PetscCall(MatDiagonalSet_Default(Y,D,is)); 205 PetscFunctionReturn(0); 206 } 207 208 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 209 { 210 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 211 PetscInt i,ishift; 212 213 PetscFunctionBegin; 214 *m = A->rmap->n; 215 if (!ia) PetscFunctionReturn(0); 216 ishift = 0; 217 if (symmetric && !A->structurally_symmetric) { 218 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja)); 219 } else if (oshift == 1) { 220 PetscInt *tia; 221 PetscInt nz = a->i[A->rmap->n]; 222 /* malloc space and add 1 to i and j indices */ 223 PetscCall(PetscMalloc1(A->rmap->n+1,&tia)); 224 for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1; 225 *ia = tia; 226 if (ja) { 227 PetscInt *tja; 228 PetscCall(PetscMalloc1(nz+1,&tja)); 229 for (i=0; i<nz; i++) tja[i] = a->j[i] + 1; 230 *ja = tja; 231 } 232 } else { 233 *ia = a->i; 234 if (ja) *ja = a->j; 235 } 236 PetscFunctionReturn(0); 237 } 238 239 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 240 { 241 PetscFunctionBegin; 242 if (!ia) PetscFunctionReturn(0); 243 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 244 PetscCall(PetscFree(*ia)); 245 if (ja) PetscCall(PetscFree(*ja)); 246 } 247 PetscFunctionReturn(0); 248 } 249 250 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 251 { 252 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 253 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 254 PetscInt nz = a->i[m],row,*jj,mr,col; 255 256 PetscFunctionBegin; 257 *nn = n; 258 if (!ia) PetscFunctionReturn(0); 259 if (symmetric) { 260 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,0,oshift,(PetscInt**)ia,(PetscInt**)ja)); 261 } else { 262 PetscCall(PetscCalloc1(n,&collengths)); 263 PetscCall(PetscMalloc1(n+1,&cia)); 264 PetscCall(PetscMalloc1(nz,&cja)); 265 jj = a->j; 266 for (i=0; i<nz; i++) { 267 collengths[jj[i]]++; 268 } 269 cia[0] = oshift; 270 for (i=0; i<n; i++) { 271 cia[i+1] = cia[i] + collengths[i]; 272 } 273 PetscCall(PetscArrayzero(collengths,n)); 274 jj = a->j; 275 for (row=0; row<m; row++) { 276 mr = a->i[row+1] - a->i[row]; 277 for (i=0; i<mr; i++) { 278 col = *jj++; 279 280 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 281 } 282 } 283 PetscCall(PetscFree(collengths)); 284 *ia = cia; *ja = cja; 285 } 286 PetscFunctionReturn(0); 287 } 288 289 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 290 { 291 PetscFunctionBegin; 292 if (!ia) PetscFunctionReturn(0); 293 294 PetscCall(PetscFree(*ia)); 295 PetscCall(PetscFree(*ja)); 296 PetscFunctionReturn(0); 297 } 298 299 /* 300 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 301 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 302 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ() 303 */ 304 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 305 { 306 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 307 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 308 PetscInt nz = a->i[m],row,mr,col,tmp; 309 PetscInt *cspidx; 310 const PetscInt *jj; 311 312 PetscFunctionBegin; 313 *nn = n; 314 if (!ia) PetscFunctionReturn(0); 315 316 PetscCall(PetscCalloc1(n,&collengths)); 317 PetscCall(PetscMalloc1(n+1,&cia)); 318 PetscCall(PetscMalloc1(nz,&cja)); 319 PetscCall(PetscMalloc1(nz,&cspidx)); 320 jj = a->j; 321 for (i=0; i<nz; i++) { 322 collengths[jj[i]]++; 323 } 324 cia[0] = oshift; 325 for (i=0; i<n; i++) { 326 cia[i+1] = cia[i] + collengths[i]; 327 } 328 PetscCall(PetscArrayzero(collengths,n)); 329 jj = a->j; 330 for (row=0; row<m; row++) { 331 mr = a->i[row+1] - a->i[row]; 332 for (i=0; i<mr; i++) { 333 col = *jj++; 334 tmp = cia[col] + collengths[col]++ - oshift; 335 cspidx[tmp] = a->i[row] + i; /* index of a->j */ 336 cja[tmp] = row + oshift; 337 } 338 } 339 PetscCall(PetscFree(collengths)); 340 *ia = cia; 341 *ja = cja; 342 *spidx = cspidx; 343 PetscFunctionReturn(0); 344 } 345 346 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 347 { 348 PetscFunctionBegin; 349 PetscCall(MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done)); 350 PetscCall(PetscFree(*spidx)); 351 PetscFunctionReturn(0); 352 } 353 354 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 355 { 356 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 357 PetscInt *ai = a->i; 358 PetscScalar *aa; 359 360 PetscFunctionBegin; 361 PetscCall(MatSeqAIJGetArray(A,&aa)); 362 PetscCall(PetscArraycpy(aa+ai[row],v,ai[row+1]-ai[row])); 363 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 364 PetscFunctionReturn(0); 365 } 366 367 /* 368 MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions 369 370 - a single row of values is set with each call 371 - no row or column indices are negative or (in error) larger than the number of rows or columns 372 - the values are always added to the matrix, not set 373 - no new locations are introduced in the nonzero structure of the matrix 374 375 This does NOT assume the global column indices are sorted 376 377 */ 378 379 #include <petsc/private/isimpl.h> 380 PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 381 { 382 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 383 PetscInt low,high,t,row,nrow,i,col,l; 384 const PetscInt *rp,*ai = a->i,*ailen = a->ilen,*aj = a->j; 385 PetscInt lastcol = -1; 386 MatScalar *ap,value,*aa; 387 const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices; 388 389 PetscFunctionBegin; 390 PetscCall(MatSeqAIJGetArray(A,&aa)); 391 row = ridx[im[0]]; 392 rp = aj + ai[row]; 393 ap = aa + ai[row]; 394 nrow = ailen[row]; 395 low = 0; 396 high = nrow; 397 for (l=0; l<n; l++) { /* loop over added columns */ 398 col = cidx[in[l]]; 399 value = v[l]; 400 401 if (col <= lastcol) low = 0; 402 else high = nrow; 403 lastcol = col; 404 while (high-low > 5) { 405 t = (low+high)/2; 406 if (rp[t] > col) high = t; 407 else low = t; 408 } 409 for (i=low; i<high; i++) { 410 if (rp[i] == col) { 411 ap[i] += value; 412 low = i + 1; 413 break; 414 } 415 } 416 } 417 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 418 return 0; 419 } 420 421 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 422 { 423 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 424 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 425 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 426 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 427 MatScalar *ap=NULL,value=0.0,*aa; 428 PetscBool ignorezeroentries = a->ignorezeroentries; 429 PetscBool roworiented = a->roworiented; 430 431 PetscFunctionBegin; 432 PetscCall(MatSeqAIJGetArray(A,&aa)); 433 for (k=0; k<m; k++) { /* loop over added rows */ 434 row = im[k]; 435 if (row < 0) continue; 436 PetscCheck(row < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->n-1); 437 rp = aj + ai[row]; 438 if (!A->structure_only) ap = aa + ai[row]; 439 rmax = imax[row]; nrow = ailen[row]; 440 low = 0; 441 high = nrow; 442 for (l=0; l<n; l++) { /* loop over added columns */ 443 if (in[l] < 0) continue; 444 PetscCheck(in[l] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->n-1); 445 col = in[l]; 446 if (v && !A->structure_only) value = roworiented ? v[l + k*n] : v[k + l*m]; 447 if (!A->structure_only && value == 0.0 && ignorezeroentries && is == ADD_VALUES && row != col) continue; 448 449 if (col <= lastcol) low = 0; 450 else high = nrow; 451 lastcol = col; 452 while (high-low > 5) { 453 t = (low+high)/2; 454 if (rp[t] > col) high = t; 455 else low = t; 456 } 457 for (i=low; i<high; i++) { 458 if (rp[i] > col) break; 459 if (rp[i] == col) { 460 if (!A->structure_only) { 461 if (is == ADD_VALUES) { 462 ap[i] += value; 463 (void)PetscLogFlops(1.0); 464 } 465 else ap[i] = value; 466 } 467 low = i + 1; 468 goto noinsert; 469 } 470 } 471 if (value == 0.0 && ignorezeroentries && row != col) goto noinsert; 472 if (nonew == 1) goto noinsert; 473 PetscCheck(nonew != -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%" PetscInt_FMT ",%" PetscInt_FMT ") in the matrix",row,col); 474 if (A->structure_only) { 475 MatSeqXAIJReallocateAIJ_structure_only(A,A->rmap->n,1,nrow,row,col,rmax,ai,aj,rp,imax,nonew,MatScalar); 476 } else { 477 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 478 } 479 N = nrow++ - 1; a->nz++; high++; 480 /* shift up all the later entries in this row */ 481 PetscCall(PetscArraymove(rp+i+1,rp+i,N-i+1)); 482 rp[i] = col; 483 if (!A->structure_only) { 484 PetscCall(PetscArraymove(ap+i+1,ap+i,N-i+1)); 485 ap[i] = value; 486 } 487 low = i + 1; 488 A->nonzerostate++; 489 noinsert:; 490 } 491 ailen[row] = nrow; 492 } 493 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 494 PetscFunctionReturn(0); 495 } 496 497 PetscErrorCode MatSetValues_SeqAIJ_SortedFullNoPreallocation(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 498 { 499 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 500 PetscInt *rp,k,row; 501 PetscInt *ai = a->i; 502 PetscInt *aj = a->j; 503 MatScalar *aa,*ap; 504 505 PetscFunctionBegin; 506 PetscCheck(!A->was_assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot call on assembled matrix."); 507 PetscCheck(m*n+a->nz <= a->maxnz,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of entries in matrix will be larger than maximum nonzeros allocated for %" PetscInt_FMT " in MatSeqAIJSetTotalPreallocation()",a->maxnz); 508 509 PetscCall(MatSeqAIJGetArray(A,&aa)); 510 for (k=0; k<m; k++) { /* loop over added rows */ 511 row = im[k]; 512 rp = aj + ai[row]; 513 ap = aa + ai[row]; 514 515 PetscCall(PetscMemcpy(rp,in,n*sizeof(PetscInt))); 516 if (!A->structure_only) { 517 if (v) { 518 PetscCall(PetscMemcpy(ap,v,n*sizeof(PetscScalar))); 519 v += n; 520 } else { 521 PetscCall(PetscMemzero(ap,n*sizeof(PetscScalar))); 522 } 523 } 524 a->ilen[row] = n; 525 a->imax[row] = n; 526 a->i[row+1] = a->i[row]+n; 527 a->nz += n; 528 } 529 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 530 PetscFunctionReturn(0); 531 } 532 533 /*@ 534 MatSeqAIJSetTotalPreallocation - Sets an upper bound on the total number of expected nonzeros in the matrix. 535 536 Input Parameters: 537 + A - the SeqAIJ matrix 538 - nztotal - bound on the number of nonzeros 539 540 Level: advanced 541 542 Notes: 543 This can be called if you will be provided the matrix row by row (from row zero) with sorted column indices for each row. 544 Simply call MatSetValues() after this call to provide the matrix entries in the usual manner. This matrix may be used 545 as always with multiple matrix assemblies. 546 547 .seealso: `MatSetOption()`, `MAT_SORTED_FULL`, `MatSetValues()`, `MatSeqAIJSetPreallocation()` 548 @*/ 549 550 PetscErrorCode MatSeqAIJSetTotalPreallocation(Mat A,PetscInt nztotal) 551 { 552 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 553 554 PetscFunctionBegin; 555 PetscCall(PetscLayoutSetUp(A->rmap)); 556 PetscCall(PetscLayoutSetUp(A->cmap)); 557 a->maxnz = nztotal; 558 if (!a->imax) { 559 PetscCall(PetscMalloc1(A->rmap->n,&a->imax)); 560 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt))); 561 } 562 if (!a->ilen) { 563 PetscCall(PetscMalloc1(A->rmap->n,&a->ilen)); 564 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt))); 565 } else { 566 PetscCall(PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt))); 567 } 568 569 /* allocate the matrix space */ 570 if (A->structure_only) { 571 PetscCall(PetscMalloc1(nztotal,&a->j)); 572 PetscCall(PetscMalloc1(A->rmap->n+1,&a->i)); 573 PetscCall(PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*sizeof(PetscInt))); 574 } else { 575 PetscCall(PetscMalloc3(nztotal,&a->a,nztotal,&a->j,A->rmap->n+1,&a->i)); 576 PetscCall(PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*(sizeof(PetscScalar)+sizeof(PetscInt)))); 577 } 578 a->i[0] = 0; 579 if (A->structure_only) { 580 a->singlemalloc = PETSC_FALSE; 581 a->free_a = PETSC_FALSE; 582 } else { 583 a->singlemalloc = PETSC_TRUE; 584 a->free_a = PETSC_TRUE; 585 } 586 a->free_ij = PETSC_TRUE; 587 A->ops->setvalues = MatSetValues_SeqAIJ_SortedFullNoPreallocation; 588 A->preallocated = PETSC_TRUE; 589 PetscFunctionReturn(0); 590 } 591 592 PetscErrorCode MatSetValues_SeqAIJ_SortedFull(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 593 { 594 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 595 PetscInt *rp,k,row; 596 PetscInt *ai = a->i,*ailen = a->ilen; 597 PetscInt *aj = a->j; 598 MatScalar *aa,*ap; 599 600 PetscFunctionBegin; 601 PetscCall(MatSeqAIJGetArray(A,&aa)); 602 for (k=0; k<m; k++) { /* loop over added rows */ 603 row = im[k]; 604 PetscCheck(n <= a->imax[row],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Preallocation for row %" PetscInt_FMT " does not match number of columns provided",n); 605 rp = aj + ai[row]; 606 ap = aa + ai[row]; 607 if (!A->was_assembled) { 608 PetscCall(PetscMemcpy(rp,in,n*sizeof(PetscInt))); 609 } 610 if (!A->structure_only) { 611 if (v) { 612 PetscCall(PetscMemcpy(ap,v,n*sizeof(PetscScalar))); 613 v += n; 614 } else { 615 PetscCall(PetscMemzero(ap,n*sizeof(PetscScalar))); 616 } 617 } 618 ailen[row] = n; 619 a->nz += n; 620 } 621 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 622 PetscFunctionReturn(0); 623 } 624 625 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 626 { 627 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 628 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 629 PetscInt *ai = a->i,*ailen = a->ilen; 630 MatScalar *ap,*aa; 631 632 PetscFunctionBegin; 633 PetscCall(MatSeqAIJGetArray(A,&aa)); 634 for (k=0; k<m; k++) { /* loop over rows */ 635 row = im[k]; 636 if (row < 0) {v += n; continue;} /* negative row */ 637 PetscCheck(row < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->n-1); 638 rp = aj + ai[row]; ap = aa + ai[row]; 639 nrow = ailen[row]; 640 for (l=0; l<n; l++) { /* loop over columns */ 641 if (in[l] < 0) {v++; continue;} /* negative column */ 642 PetscCheck(in[l] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->n-1); 643 col = in[l]; 644 high = nrow; low = 0; /* assume unsorted */ 645 while (high-low > 5) { 646 t = (low+high)/2; 647 if (rp[t] > col) high = t; 648 else low = t; 649 } 650 for (i=low; i<high; i++) { 651 if (rp[i] > col) break; 652 if (rp[i] == col) { 653 *v++ = ap[i]; 654 goto finished; 655 } 656 } 657 *v++ = 0.0; 658 finished:; 659 } 660 } 661 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 662 PetscFunctionReturn(0); 663 } 664 665 PetscErrorCode MatView_SeqAIJ_Binary(Mat mat,PetscViewer viewer) 666 { 667 Mat_SeqAIJ *A = (Mat_SeqAIJ*)mat->data; 668 const PetscScalar *av; 669 PetscInt header[4],M,N,m,nz,i; 670 PetscInt *rowlens; 671 672 PetscFunctionBegin; 673 PetscCall(PetscViewerSetUp(viewer)); 674 675 M = mat->rmap->N; 676 N = mat->cmap->N; 677 m = mat->rmap->n; 678 nz = A->nz; 679 680 /* write matrix header */ 681 header[0] = MAT_FILE_CLASSID; 682 header[1] = M; header[2] = N; header[3] = nz; 683 PetscCall(PetscViewerBinaryWrite(viewer,header,4,PETSC_INT)); 684 685 /* fill in and store row lengths */ 686 PetscCall(PetscMalloc1(m,&rowlens)); 687 for (i=0; i<m; i++) rowlens[i] = A->i[i+1] - A->i[i]; 688 PetscCall(PetscViewerBinaryWrite(viewer,rowlens,m,PETSC_INT)); 689 PetscCall(PetscFree(rowlens)); 690 /* store column indices */ 691 PetscCall(PetscViewerBinaryWrite(viewer,A->j,nz,PETSC_INT)); 692 /* store nonzero values */ 693 PetscCall(MatSeqAIJGetArrayRead(mat,&av)); 694 PetscCall(PetscViewerBinaryWrite(viewer,av,nz,PETSC_SCALAR)); 695 PetscCall(MatSeqAIJRestoreArrayRead(mat,&av)); 696 697 /* write block size option to the viewer's .info file */ 698 PetscCall(MatView_Binary_BlockSizes(mat,viewer)); 699 PetscFunctionReturn(0); 700 } 701 702 static PetscErrorCode MatView_SeqAIJ_ASCII_structonly(Mat A,PetscViewer viewer) 703 { 704 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 705 PetscInt i,k,m=A->rmap->N; 706 707 PetscFunctionBegin; 708 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 709 for (i=0; i<m; i++) { 710 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 711 for (k=a->i[i]; k<a->i[i+1]; k++) { 712 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ") ",a->j[k])); 713 } 714 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 715 } 716 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 717 PetscFunctionReturn(0); 718 } 719 720 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 721 722 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 723 { 724 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 725 const PetscScalar *av; 726 PetscInt i,j,m = A->rmap->n; 727 const char *name; 728 PetscViewerFormat format; 729 730 PetscFunctionBegin; 731 if (A->structure_only) { 732 PetscCall(MatView_SeqAIJ_ASCII_structonly(A,viewer)); 733 PetscFunctionReturn(0); 734 } 735 736 PetscCall(PetscViewerGetFormat(viewer,&format)); 737 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 738 739 /* trigger copy to CPU if needed */ 740 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 741 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 742 if (format == PETSC_VIEWER_ASCII_MATLAB) { 743 PetscInt nofinalvalue = 0; 744 if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) { 745 /* Need a dummy value to ensure the dimension of the matrix. */ 746 nofinalvalue = 1; 747 } 748 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 749 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n",m,A->cmap->n)); 750 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %" PetscInt_FMT " \n",a->nz)); 751 #if defined(PETSC_USE_COMPLEX) 752 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",4);\n",a->nz+nofinalvalue)); 753 #else 754 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",a->nz+nofinalvalue)); 755 #endif 756 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n")); 757 758 for (i=0; i<m; i++) { 759 for (j=a->i[i]; j<a->i[i+1]; j++) { 760 #if defined(PETSC_USE_COMPLEX) 761 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",i+1,a->j[j]+1,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 762 #else 763 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,a->j[j]+1,(double)a->a[j])); 764 #endif 765 } 766 } 767 if (nofinalvalue) { 768 #if defined(PETSC_USE_COMPLEX) 769 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",m,A->cmap->n,0.,0.)); 770 #else 771 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",m,A->cmap->n,0.0)); 772 #endif 773 } 774 PetscCall(PetscObjectGetName((PetscObject)A,&name)); 775 PetscCall(PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name)); 776 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 777 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 778 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 779 for (i=0; i<m; i++) { 780 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 781 for (j=a->i[i]; j<a->i[i+1]; j++) { 782 #if defined(PETSC_USE_COMPLEX) 783 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 784 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 785 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 786 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]))); 787 } else if (PetscRealPart(a->a[j]) != 0.0) { 788 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 789 } 790 #else 791 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 792 #endif 793 } 794 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 795 } 796 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 797 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 798 PetscInt nzd=0,fshift=1,*sptr; 799 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 800 PetscCall(PetscMalloc1(m+1,&sptr)); 801 for (i=0; i<m; i++) { 802 sptr[i] = nzd+1; 803 for (j=a->i[i]; j<a->i[i+1]; j++) { 804 if (a->j[j] >= i) { 805 #if defined(PETSC_USE_COMPLEX) 806 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 807 #else 808 if (a->a[j] != 0.0) nzd++; 809 #endif 810 } 811 } 812 } 813 sptr[m] = nzd+1; 814 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n\n",m,nzd)); 815 for (i=0; i<m+1; i+=6) { 816 if (i+4<m) { 817 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5])); 818 } else if (i+3<m) { 819 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4])); 820 } else if (i+2<m) { 821 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3])); 822 } else if (i+1<m) { 823 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2])); 824 } else if (i<m) { 825 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1])); 826 } else { 827 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT "\n",sptr[i])); 828 } 829 } 830 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 831 PetscCall(PetscFree(sptr)); 832 for (i=0; i<m; i++) { 833 for (j=a->i[i]; j<a->i[i+1]; j++) { 834 if (a->j[j] >= i) PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " ",a->j[j]+fshift)); 835 } 836 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 837 } 838 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 839 for (i=0; i<m; i++) { 840 for (j=a->i[i]; j<a->i[i+1]; j++) { 841 if (a->j[j] >= i) { 842 #if defined(PETSC_USE_COMPLEX) 843 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 844 PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 845 } 846 #else 847 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)a->a[j])); 848 #endif 849 } 850 } 851 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 852 } 853 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 854 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 855 PetscInt cnt = 0,jcnt; 856 PetscScalar value; 857 #if defined(PETSC_USE_COMPLEX) 858 PetscBool realonly = PETSC_TRUE; 859 860 for (i=0; i<a->i[m]; i++) { 861 if (PetscImaginaryPart(a->a[i]) != 0.0) { 862 realonly = PETSC_FALSE; 863 break; 864 } 865 } 866 #endif 867 868 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 869 for (i=0; i<m; i++) { 870 jcnt = 0; 871 for (j=0; j<A->cmap->n; j++) { 872 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 873 value = a->a[cnt++]; 874 jcnt++; 875 } else { 876 value = 0.0; 877 } 878 #if defined(PETSC_USE_COMPLEX) 879 if (realonly) { 880 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e ",(double)PetscRealPart(value))); 881 } else { 882 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",(double)PetscRealPart(value),(double)PetscImaginaryPart(value))); 883 } 884 #else 885 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e ",(double)value)); 886 #endif 887 } 888 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 889 } 890 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 891 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 892 PetscInt fshift=1; 893 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 894 #if defined(PETSC_USE_COMPLEX) 895 PetscCall(PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate complex general\n")); 896 #else 897 PetscCall(PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate real general\n")); 898 #endif 899 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", m, A->cmap->n, a->nz)); 900 for (i=0; i<m; i++) { 901 for (j=a->i[i]; j<a->i[i+1]; j++) { 902 #if defined(PETSC_USE_COMPLEX) 903 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g %g\n", i+fshift,a->j[j]+fshift,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 904 #else 905 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j])); 906 #endif 907 } 908 } 909 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 910 } else { 911 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 912 if (A->factortype) { 913 for (i=0; i<m; i++) { 914 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 915 /* L part */ 916 for (j=a->i[i]; j<a->i[i+1]; j++) { 917 #if defined(PETSC_USE_COMPLEX) 918 if (PetscImaginaryPart(a->a[j]) > 0.0) { 919 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 920 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 921 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])))); 922 } else { 923 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 924 } 925 #else 926 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 927 #endif 928 } 929 /* diagonal */ 930 j = a->diag[i]; 931 #if defined(PETSC_USE_COMPLEX) 932 if (PetscImaginaryPart(a->a[j]) > 0.0) { 933 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)PetscImaginaryPart(1.0/a->a[j]))); 934 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 935 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)(-PetscImaginaryPart(1.0/a->a[j])))); 936 } else { 937 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(1.0/a->a[j]))); 938 } 939 #else 940 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)(1.0/a->a[j]))); 941 #endif 942 943 /* U part */ 944 for (j=a->diag[i+1]+1; j<a->diag[i]; j++) { 945 #if defined(PETSC_USE_COMPLEX) 946 if (PetscImaginaryPart(a->a[j]) > 0.0) { 947 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 948 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 949 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])))); 950 } else { 951 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 952 } 953 #else 954 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 955 #endif 956 } 957 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 958 } 959 } else { 960 for (i=0; i<m; i++) { 961 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 962 for (j=a->i[i]; j<a->i[i+1]; j++) { 963 #if defined(PETSC_USE_COMPLEX) 964 if (PetscImaginaryPart(a->a[j]) > 0.0) { 965 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 966 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 967 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]))); 968 } else { 969 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 970 } 971 #else 972 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 973 #endif 974 } 975 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 976 } 977 } 978 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 979 } 980 PetscCall(PetscViewerFlush(viewer)); 981 PetscFunctionReturn(0); 982 } 983 984 #include <petscdraw.h> 985 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 986 { 987 Mat A = (Mat) Aa; 988 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 989 PetscInt i,j,m = A->rmap->n; 990 int color; 991 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 992 PetscViewer viewer; 993 PetscViewerFormat format; 994 const PetscScalar *aa; 995 996 PetscFunctionBegin; 997 PetscCall(PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer)); 998 PetscCall(PetscViewerGetFormat(viewer,&format)); 999 PetscCall(PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr)); 1000 1001 /* loop over matrix elements drawing boxes */ 1002 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1003 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1004 PetscDrawCollectiveBegin(draw); 1005 /* Blue for negative, Cyan for zero and Red for positive */ 1006 color = PETSC_DRAW_BLUE; 1007 for (i=0; i<m; i++) { 1008 y_l = m - i - 1.0; y_r = y_l + 1.0; 1009 for (j=a->i[i]; j<a->i[i+1]; j++) { 1010 x_l = a->j[j]; x_r = x_l + 1.0; 1011 if (PetscRealPart(aa[j]) >= 0.) continue; 1012 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1013 } 1014 } 1015 color = PETSC_DRAW_CYAN; 1016 for (i=0; i<m; i++) { 1017 y_l = m - i - 1.0; y_r = y_l + 1.0; 1018 for (j=a->i[i]; j<a->i[i+1]; j++) { 1019 x_l = a->j[j]; x_r = x_l + 1.0; 1020 if (aa[j] != 0.) continue; 1021 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1022 } 1023 } 1024 color = PETSC_DRAW_RED; 1025 for (i=0; i<m; i++) { 1026 y_l = m - i - 1.0; y_r = y_l + 1.0; 1027 for (j=a->i[i]; j<a->i[i+1]; j++) { 1028 x_l = a->j[j]; x_r = x_l + 1.0; 1029 if (PetscRealPart(aa[j]) <= 0.) continue; 1030 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1031 } 1032 } 1033 PetscDrawCollectiveEnd(draw); 1034 } else { 1035 /* use contour shading to indicate magnitude of values */ 1036 /* first determine max of all nonzero values */ 1037 PetscReal minv = 0.0, maxv = 0.0; 1038 PetscInt nz = a->nz, count = 0; 1039 PetscDraw popup; 1040 1041 for (i=0; i<nz; i++) { 1042 if (PetscAbsScalar(aa[i]) > maxv) maxv = PetscAbsScalar(aa[i]); 1043 } 1044 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1045 PetscCall(PetscDrawGetPopup(draw,&popup)); 1046 PetscCall(PetscDrawScalePopup(popup,minv,maxv)); 1047 1048 PetscDrawCollectiveBegin(draw); 1049 for (i=0; i<m; i++) { 1050 y_l = m - i - 1.0; 1051 y_r = y_l + 1.0; 1052 for (j=a->i[i]; j<a->i[i+1]; j++) { 1053 x_l = a->j[j]; 1054 x_r = x_l + 1.0; 1055 color = PetscDrawRealToColor(PetscAbsScalar(aa[count]),minv,maxv); 1056 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1057 count++; 1058 } 1059 } 1060 PetscDrawCollectiveEnd(draw); 1061 } 1062 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1063 PetscFunctionReturn(0); 1064 } 1065 1066 #include <petscdraw.h> 1067 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 1068 { 1069 PetscDraw draw; 1070 PetscReal xr,yr,xl,yl,h,w; 1071 PetscBool isnull; 1072 1073 PetscFunctionBegin; 1074 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 1075 PetscCall(PetscDrawIsNull(draw,&isnull)); 1076 if (isnull) PetscFunctionReturn(0); 1077 1078 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1079 xr += w; yr += h; xl = -w; yl = -h; 1080 PetscCall(PetscDrawSetCoordinates(draw,xl,yl,xr,yr)); 1081 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer)); 1082 PetscCall(PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A)); 1083 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL)); 1084 PetscCall(PetscDrawSave(draw)); 1085 PetscFunctionReturn(0); 1086 } 1087 1088 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 1089 { 1090 PetscBool iascii,isbinary,isdraw; 1091 1092 PetscFunctionBegin; 1093 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 1094 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 1095 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 1096 if (iascii) PetscCall(MatView_SeqAIJ_ASCII(A,viewer)); 1097 else if (isbinary) PetscCall(MatView_SeqAIJ_Binary(A,viewer)); 1098 else if (isdraw) PetscCall(MatView_SeqAIJ_Draw(A,viewer)); 1099 PetscCall(MatView_SeqAIJ_Inode(A,viewer)); 1100 PetscFunctionReturn(0); 1101 } 1102 1103 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 1104 { 1105 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1106 PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax; 1107 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 1108 MatScalar *aa = a->a,*ap; 1109 PetscReal ratio = 0.6; 1110 1111 PetscFunctionBegin; 1112 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 1113 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1114 if (A->was_assembled && A->ass_nonzerostate == A->nonzerostate) { 1115 /* we need to respect users asking to use or not the inodes routine in between matrix assemblies */ 1116 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A,mode)); 1117 PetscFunctionReturn(0); 1118 } 1119 1120 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 1121 for (i=1; i<m; i++) { 1122 /* move each row back by the amount of empty slots (fshift) before it*/ 1123 fshift += imax[i-1] - ailen[i-1]; 1124 rmax = PetscMax(rmax,ailen[i]); 1125 if (fshift) { 1126 ip = aj + ai[i]; 1127 ap = aa + ai[i]; 1128 N = ailen[i]; 1129 PetscCall(PetscArraymove(ip-fshift,ip,N)); 1130 if (!A->structure_only) { 1131 PetscCall(PetscArraymove(ap-fshift,ap,N)); 1132 } 1133 } 1134 ai[i] = ai[i-1] + ailen[i-1]; 1135 } 1136 if (m) { 1137 fshift += imax[m-1] - ailen[m-1]; 1138 ai[m] = ai[m-1] + ailen[m-1]; 1139 } 1140 1141 /* reset ilen and imax for each row */ 1142 a->nonzerorowcnt = 0; 1143 if (A->structure_only) { 1144 PetscCall(PetscFree(a->imax)); 1145 PetscCall(PetscFree(a->ilen)); 1146 } else { /* !A->structure_only */ 1147 for (i=0; i<m; i++) { 1148 ailen[i] = imax[i] = ai[i+1] - ai[i]; 1149 a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0); 1150 } 1151 } 1152 a->nz = ai[m]; 1153 PetscCheck(!fshift || a->nounused != -1,PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %" PetscInt_FMT " X %" PetscInt_FMT ", %" PetscInt_FMT " unneeded", m, A->cmap->n, fshift); 1154 1155 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1156 PetscCall(PetscInfo(A,"Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded,%" PetscInt_FMT " used\n",m,A->cmap->n,fshift,a->nz)); 1157 PetscCall(PetscInfo(A,"Number of mallocs during MatSetValues() is %" PetscInt_FMT "\n",a->reallocs)); 1158 PetscCall(PetscInfo(A,"Maximum nonzeros in any row is %" PetscInt_FMT "\n",rmax)); 1159 1160 A->info.mallocs += a->reallocs; 1161 a->reallocs = 0; 1162 A->info.nz_unneeded = (PetscReal)fshift; 1163 a->rmax = rmax; 1164 1165 if (!A->structure_only) { 1166 PetscCall(MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio)); 1167 } 1168 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A,mode)); 1169 PetscFunctionReturn(0); 1170 } 1171 1172 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 1173 { 1174 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1175 PetscInt i,nz = a->nz; 1176 MatScalar *aa; 1177 1178 PetscFunctionBegin; 1179 PetscCall(MatSeqAIJGetArray(A,&aa)); 1180 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1181 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 1182 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1183 PetscFunctionReturn(0); 1184 } 1185 1186 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 1187 { 1188 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1189 PetscInt i,nz = a->nz; 1190 MatScalar *aa; 1191 1192 PetscFunctionBegin; 1193 PetscCall(MatSeqAIJGetArray(A,&aa)); 1194 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1195 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 1196 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1197 PetscFunctionReturn(0); 1198 } 1199 1200 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1201 { 1202 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1203 MatScalar *aa; 1204 1205 PetscFunctionBegin; 1206 PetscCall(MatSeqAIJGetArrayWrite(A,&aa)); 1207 PetscCall(PetscArrayzero(aa,a->i[A->rmap->n])); 1208 PetscCall(MatSeqAIJRestoreArrayWrite(A,&aa)); 1209 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1210 PetscFunctionReturn(0); 1211 } 1212 1213 PETSC_INTERN PetscErrorCode MatResetPreallocationCOO_SeqAIJ(Mat A) 1214 { 1215 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1216 1217 PetscFunctionBegin; 1218 PetscCall(PetscFree(a->perm)); 1219 PetscCall(PetscFree(a->jmap)); 1220 PetscFunctionReturn(0); 1221 } 1222 1223 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1224 { 1225 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1226 1227 PetscFunctionBegin; 1228 #if defined(PETSC_USE_LOG) 1229 PetscLogObjectState((PetscObject)A,"Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT,A->rmap->n,A->cmap->n,a->nz); 1230 #endif 1231 PetscCall(MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i)); 1232 PetscCall(MatResetPreallocationCOO_SeqAIJ(A)); 1233 PetscCall(ISDestroy(&a->row)); 1234 PetscCall(ISDestroy(&a->col)); 1235 PetscCall(PetscFree(a->diag)); 1236 PetscCall(PetscFree(a->ibdiag)); 1237 PetscCall(PetscFree(a->imax)); 1238 PetscCall(PetscFree(a->ilen)); 1239 PetscCall(PetscFree(a->ipre)); 1240 PetscCall(PetscFree3(a->idiag,a->mdiag,a->ssor_work)); 1241 PetscCall(PetscFree(a->solve_work)); 1242 PetscCall(ISDestroy(&a->icol)); 1243 PetscCall(PetscFree(a->saved_values)); 1244 PetscCall(PetscFree2(a->compressedrow.i,a->compressedrow.rindex)); 1245 PetscCall(MatDestroy_SeqAIJ_Inode(A)); 1246 PetscCall(PetscFree(A->data)); 1247 1248 /* MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted may allocate this. 1249 That function is so heavily used (sometimes in an hidden way through multnumeric function pointers) 1250 that is hard to properly add this data to the MatProduct data. We free it here to avoid 1251 users reusing the matrix object with different data to incur in obscure segmentation faults 1252 due to different matrix sizes */ 1253 PetscCall(PetscObjectCompose((PetscObject)A,"__PETSc__ab_dense",NULL)); 1254 1255 PetscCall(PetscObjectChangeTypeName((PetscObject)A,NULL)); 1256 PetscCall(PetscObjectComposeFunction((PetscObject)A,"PetscMatlabEnginePut_C",NULL)); 1257 PetscCall(PetscObjectComposeFunction((PetscObject)A,"PetscMatlabEngineGet_C",NULL)); 1258 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL)); 1259 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL)); 1260 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL)); 1261 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL)); 1262 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL)); 1263 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL)); 1264 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijsell_C",NULL)); 1265 #if defined(PETSC_HAVE_MKL_SPARSE) 1266 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijmkl_C",NULL)); 1267 #endif 1268 #if defined(PETSC_HAVE_CUDA) 1269 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcusparse_C",NULL)); 1270 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",NULL)); 1271 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",NULL)); 1272 #endif 1273 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 1274 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijkokkos_C",NULL)); 1275 #endif 1276 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcrl_C",NULL)); 1277 #if defined(PETSC_HAVE_ELEMENTAL) 1278 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_elemental_C",NULL)); 1279 #endif 1280 #if defined(PETSC_HAVE_SCALAPACK) 1281 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_scalapack_C",NULL)); 1282 #endif 1283 #if defined(PETSC_HAVE_HYPRE) 1284 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_hypre_C",NULL)); 1285 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",NULL)); 1286 #endif 1287 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqdense_C",NULL)); 1288 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsell_C",NULL)); 1289 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_is_C",NULL)); 1290 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL)); 1291 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatIsHermitianTranspose_C",NULL)); 1292 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL)); 1293 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatResetPreallocation_C",NULL)); 1294 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL)); 1295 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL)); 1296 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_is_seqaij_C",NULL)); 1297 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqdense_seqaij_C",NULL)); 1298 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaij_C",NULL)); 1299 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJKron_C",NULL)); 1300 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL)); 1301 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",NULL)); 1302 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL)); 1303 /* these calls do not belong here: the subclasses Duplicate/Destroy are wrong */ 1304 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaijsell_seqaij_C",NULL)); 1305 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaijperm_seqaij_C",NULL)); 1306 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijviennacl_C",NULL)); 1307 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaijviennacl_seqdense_C",NULL)); 1308 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaijviennacl_seqaij_C",NULL)); 1309 PetscFunctionReturn(0); 1310 } 1311 1312 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1313 { 1314 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1315 1316 PetscFunctionBegin; 1317 switch (op) { 1318 case MAT_ROW_ORIENTED: 1319 a->roworiented = flg; 1320 break; 1321 case MAT_KEEP_NONZERO_PATTERN: 1322 a->keepnonzeropattern = flg; 1323 break; 1324 case MAT_NEW_NONZERO_LOCATIONS: 1325 a->nonew = (flg ? 0 : 1); 1326 break; 1327 case MAT_NEW_NONZERO_LOCATION_ERR: 1328 a->nonew = (flg ? -1 : 0); 1329 break; 1330 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1331 a->nonew = (flg ? -2 : 0); 1332 break; 1333 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1334 a->nounused = (flg ? -1 : 0); 1335 break; 1336 case MAT_IGNORE_ZERO_ENTRIES: 1337 a->ignorezeroentries = flg; 1338 break; 1339 case MAT_SPD: 1340 case MAT_SYMMETRIC: 1341 case MAT_STRUCTURALLY_SYMMETRIC: 1342 case MAT_HERMITIAN: 1343 case MAT_SYMMETRY_ETERNAL: 1344 case MAT_STRUCTURE_ONLY: 1345 /* These options are handled directly by MatSetOption() */ 1346 break; 1347 case MAT_FORCE_DIAGONAL_ENTRIES: 1348 case MAT_IGNORE_OFF_PROC_ENTRIES: 1349 case MAT_USE_HASH_TABLE: 1350 PetscCall(PetscInfo(A,"Option %s ignored\n",MatOptions[op])); 1351 break; 1352 case MAT_USE_INODES: 1353 PetscCall(MatSetOption_SeqAIJ_Inode(A,MAT_USE_INODES,flg)); 1354 break; 1355 case MAT_SUBMAT_SINGLEIS: 1356 A->submat_singleis = flg; 1357 break; 1358 case MAT_SORTED_FULL: 1359 if (flg) A->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 1360 else A->ops->setvalues = MatSetValues_SeqAIJ; 1361 break; 1362 case MAT_FORM_EXPLICIT_TRANSPOSE: 1363 A->form_explicit_transpose = flg; 1364 break; 1365 default: 1366 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1367 } 1368 PetscFunctionReturn(0); 1369 } 1370 1371 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1372 { 1373 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1374 PetscInt i,j,n,*ai=a->i,*aj=a->j; 1375 PetscScalar *x; 1376 const PetscScalar *aa; 1377 1378 PetscFunctionBegin; 1379 PetscCall(VecGetLocalSize(v,&n)); 1380 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1381 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1382 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1383 PetscInt *diag=a->diag; 1384 PetscCall(VecGetArrayWrite(v,&x)); 1385 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1386 PetscCall(VecRestoreArrayWrite(v,&x)); 1387 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1388 PetscFunctionReturn(0); 1389 } 1390 1391 PetscCall(VecGetArrayWrite(v,&x)); 1392 for (i=0; i<n; i++) { 1393 x[i] = 0.0; 1394 for (j=ai[i]; j<ai[i+1]; j++) { 1395 if (aj[j] == i) { 1396 x[i] = aa[j]; 1397 break; 1398 } 1399 } 1400 } 1401 PetscCall(VecRestoreArrayWrite(v,&x)); 1402 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1403 PetscFunctionReturn(0); 1404 } 1405 1406 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1407 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1408 { 1409 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1410 const MatScalar *aa; 1411 PetscScalar *y; 1412 const PetscScalar *x; 1413 PetscInt m = A->rmap->n; 1414 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1415 const MatScalar *v; 1416 PetscScalar alpha; 1417 PetscInt n,i,j; 1418 const PetscInt *idx,*ii,*ridx=NULL; 1419 Mat_CompressedRow cprow = a->compressedrow; 1420 PetscBool usecprow = cprow.use; 1421 #endif 1422 1423 PetscFunctionBegin; 1424 if (zz != yy) PetscCall(VecCopy(zz,yy)); 1425 PetscCall(VecGetArrayRead(xx,&x)); 1426 PetscCall(VecGetArray(yy,&y)); 1427 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1428 1429 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1430 fortranmulttransposeaddaij_(&m,x,a->i,a->j,aa,y); 1431 #else 1432 if (usecprow) { 1433 m = cprow.nrows; 1434 ii = cprow.i; 1435 ridx = cprow.rindex; 1436 } else { 1437 ii = a->i; 1438 } 1439 for (i=0; i<m; i++) { 1440 idx = a->j + ii[i]; 1441 v = aa + ii[i]; 1442 n = ii[i+1] - ii[i]; 1443 if (usecprow) { 1444 alpha = x[ridx[i]]; 1445 } else { 1446 alpha = x[i]; 1447 } 1448 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1449 } 1450 #endif 1451 PetscCall(PetscLogFlops(2.0*a->nz)); 1452 PetscCall(VecRestoreArrayRead(xx,&x)); 1453 PetscCall(VecRestoreArray(yy,&y)); 1454 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1455 PetscFunctionReturn(0); 1456 } 1457 1458 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1459 { 1460 PetscFunctionBegin; 1461 PetscCall(VecSet(yy,0.0)); 1462 PetscCall(MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy)); 1463 PetscFunctionReturn(0); 1464 } 1465 1466 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1467 1468 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1469 { 1470 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1471 PetscScalar *y; 1472 const PetscScalar *x; 1473 const MatScalar *aa,*a_a; 1474 PetscInt m=A->rmap->n; 1475 const PetscInt *aj,*ii,*ridx=NULL; 1476 PetscInt n,i; 1477 PetscScalar sum; 1478 PetscBool usecprow=a->compressedrow.use; 1479 1480 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1481 #pragma disjoint(*x,*y,*aa) 1482 #endif 1483 1484 PetscFunctionBegin; 1485 if (a->inode.use && a->inode.checked) { 1486 PetscCall(MatMult_SeqAIJ_Inode(A,xx,yy)); 1487 PetscFunctionReturn(0); 1488 } 1489 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1490 PetscCall(VecGetArrayRead(xx,&x)); 1491 PetscCall(VecGetArray(yy,&y)); 1492 ii = a->i; 1493 if (usecprow) { /* use compressed row format */ 1494 PetscCall(PetscArrayzero(y,m)); 1495 m = a->compressedrow.nrows; 1496 ii = a->compressedrow.i; 1497 ridx = a->compressedrow.rindex; 1498 for (i=0; i<m; i++) { 1499 n = ii[i+1] - ii[i]; 1500 aj = a->j + ii[i]; 1501 aa = a_a + ii[i]; 1502 sum = 0.0; 1503 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1504 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1505 y[*ridx++] = sum; 1506 } 1507 } else { /* do not use compressed row format */ 1508 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1509 aj = a->j; 1510 aa = a_a; 1511 fortranmultaij_(&m,x,ii,aj,aa,y); 1512 #else 1513 for (i=0; i<m; i++) { 1514 n = ii[i+1] - ii[i]; 1515 aj = a->j + ii[i]; 1516 aa = a_a + ii[i]; 1517 sum = 0.0; 1518 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1519 y[i] = sum; 1520 } 1521 #endif 1522 } 1523 PetscCall(PetscLogFlops(2.0*a->nz - a->nonzerorowcnt)); 1524 PetscCall(VecRestoreArrayRead(xx,&x)); 1525 PetscCall(VecRestoreArray(yy,&y)); 1526 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1527 PetscFunctionReturn(0); 1528 } 1529 1530 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1531 { 1532 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1533 PetscScalar *y; 1534 const PetscScalar *x; 1535 const MatScalar *aa,*a_a; 1536 PetscInt m=A->rmap->n; 1537 const PetscInt *aj,*ii,*ridx=NULL; 1538 PetscInt n,i,nonzerorow=0; 1539 PetscScalar sum; 1540 PetscBool usecprow=a->compressedrow.use; 1541 1542 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1543 #pragma disjoint(*x,*y,*aa) 1544 #endif 1545 1546 PetscFunctionBegin; 1547 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1548 PetscCall(VecGetArrayRead(xx,&x)); 1549 PetscCall(VecGetArray(yy,&y)); 1550 if (usecprow) { /* use compressed row format */ 1551 m = a->compressedrow.nrows; 1552 ii = a->compressedrow.i; 1553 ridx = a->compressedrow.rindex; 1554 for (i=0; i<m; i++) { 1555 n = ii[i+1] - ii[i]; 1556 aj = a->j + ii[i]; 1557 aa = a_a + ii[i]; 1558 sum = 0.0; 1559 nonzerorow += (n>0); 1560 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1561 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1562 y[*ridx++] = sum; 1563 } 1564 } else { /* do not use compressed row format */ 1565 ii = a->i; 1566 for (i=0; i<m; i++) { 1567 n = ii[i+1] - ii[i]; 1568 aj = a->j + ii[i]; 1569 aa = a_a + ii[i]; 1570 sum = 0.0; 1571 nonzerorow += (n>0); 1572 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1573 y[i] = sum; 1574 } 1575 } 1576 PetscCall(PetscLogFlops(2.0*a->nz - nonzerorow)); 1577 PetscCall(VecRestoreArrayRead(xx,&x)); 1578 PetscCall(VecRestoreArray(yy,&y)); 1579 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1580 PetscFunctionReturn(0); 1581 } 1582 1583 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1584 { 1585 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1586 PetscScalar *y,*z; 1587 const PetscScalar *x; 1588 const MatScalar *aa,*a_a; 1589 PetscInt m = A->rmap->n,*aj,*ii; 1590 PetscInt n,i,*ridx=NULL; 1591 PetscScalar sum; 1592 PetscBool usecprow=a->compressedrow.use; 1593 1594 PetscFunctionBegin; 1595 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1596 PetscCall(VecGetArrayRead(xx,&x)); 1597 PetscCall(VecGetArrayPair(yy,zz,&y,&z)); 1598 if (usecprow) { /* use compressed row format */ 1599 if (zz != yy) { 1600 PetscCall(PetscArraycpy(z,y,m)); 1601 } 1602 m = a->compressedrow.nrows; 1603 ii = a->compressedrow.i; 1604 ridx = a->compressedrow.rindex; 1605 for (i=0; i<m; i++) { 1606 n = ii[i+1] - ii[i]; 1607 aj = a->j + ii[i]; 1608 aa = a_a + ii[i]; 1609 sum = y[*ridx]; 1610 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1611 z[*ridx++] = sum; 1612 } 1613 } else { /* do not use compressed row format */ 1614 ii = a->i; 1615 for (i=0; i<m; i++) { 1616 n = ii[i+1] - ii[i]; 1617 aj = a->j + ii[i]; 1618 aa = a_a + ii[i]; 1619 sum = y[i]; 1620 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1621 z[i] = sum; 1622 } 1623 } 1624 PetscCall(PetscLogFlops(2.0*a->nz)); 1625 PetscCall(VecRestoreArrayRead(xx,&x)); 1626 PetscCall(VecRestoreArrayPair(yy,zz,&y,&z)); 1627 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1628 PetscFunctionReturn(0); 1629 } 1630 1631 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1632 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1633 { 1634 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1635 PetscScalar *y,*z; 1636 const PetscScalar *x; 1637 const MatScalar *aa,*a_a; 1638 const PetscInt *aj,*ii,*ridx=NULL; 1639 PetscInt m = A->rmap->n,n,i; 1640 PetscScalar sum; 1641 PetscBool usecprow=a->compressedrow.use; 1642 1643 PetscFunctionBegin; 1644 if (a->inode.use && a->inode.checked) { 1645 PetscCall(MatMultAdd_SeqAIJ_Inode(A,xx,yy,zz)); 1646 PetscFunctionReturn(0); 1647 } 1648 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1649 PetscCall(VecGetArrayRead(xx,&x)); 1650 PetscCall(VecGetArrayPair(yy,zz,&y,&z)); 1651 if (usecprow) { /* use compressed row format */ 1652 if (zz != yy) { 1653 PetscCall(PetscArraycpy(z,y,m)); 1654 } 1655 m = a->compressedrow.nrows; 1656 ii = a->compressedrow.i; 1657 ridx = a->compressedrow.rindex; 1658 for (i=0; i<m; i++) { 1659 n = ii[i+1] - ii[i]; 1660 aj = a->j + ii[i]; 1661 aa = a_a + ii[i]; 1662 sum = y[*ridx]; 1663 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1664 z[*ridx++] = sum; 1665 } 1666 } else { /* do not use compressed row format */ 1667 ii = a->i; 1668 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1669 aj = a->j; 1670 aa = a_a; 1671 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1672 #else 1673 for (i=0; i<m; i++) { 1674 n = ii[i+1] - ii[i]; 1675 aj = a->j + ii[i]; 1676 aa = a_a + ii[i]; 1677 sum = y[i]; 1678 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1679 z[i] = sum; 1680 } 1681 #endif 1682 } 1683 PetscCall(PetscLogFlops(2.0*a->nz)); 1684 PetscCall(VecRestoreArrayRead(xx,&x)); 1685 PetscCall(VecRestoreArrayPair(yy,zz,&y,&z)); 1686 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1687 PetscFunctionReturn(0); 1688 } 1689 1690 /* 1691 Adds diagonal pointers to sparse matrix structure. 1692 */ 1693 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1694 { 1695 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1696 PetscInt i,j,m = A->rmap->n; 1697 PetscBool alreadySet = PETSC_TRUE; 1698 1699 PetscFunctionBegin; 1700 if (!a->diag) { 1701 PetscCall(PetscMalloc1(m,&a->diag)); 1702 PetscCall(PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt))); 1703 alreadySet = PETSC_FALSE; 1704 } 1705 for (i=0; i<A->rmap->n; i++) { 1706 /* If A's diagonal is already correctly set, this fast track enables cheap and repeated MatMarkDiagonal_SeqAIJ() calls */ 1707 if (alreadySet) { 1708 PetscInt pos = a->diag[i]; 1709 if (pos >= a->i[i] && pos < a->i[i+1] && a->j[pos] == i) continue; 1710 } 1711 1712 a->diag[i] = a->i[i+1]; 1713 for (j=a->i[i]; j<a->i[i+1]; j++) { 1714 if (a->j[j] == i) { 1715 a->diag[i] = j; 1716 break; 1717 } 1718 } 1719 } 1720 PetscFunctionReturn(0); 1721 } 1722 1723 PetscErrorCode MatShift_SeqAIJ(Mat A,PetscScalar v) 1724 { 1725 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1726 const PetscInt *diag = (const PetscInt*)a->diag; 1727 const PetscInt *ii = (const PetscInt*) a->i; 1728 PetscInt i,*mdiag = NULL; 1729 PetscInt cnt = 0; /* how many diagonals are missing */ 1730 1731 PetscFunctionBegin; 1732 if (!A->preallocated || !a->nz) { 1733 PetscCall(MatSeqAIJSetPreallocation(A,1,NULL)); 1734 PetscCall(MatShift_Basic(A,v)); 1735 PetscFunctionReturn(0); 1736 } 1737 1738 if (a->diagonaldense) { 1739 cnt = 0; 1740 } else { 1741 PetscCall(PetscCalloc1(A->rmap->n,&mdiag)); 1742 for (i=0; i<A->rmap->n; i++) { 1743 if (i < A->cmap->n && diag[i] >= ii[i+1]) { /* 'out of range' rows never have diagonals */ 1744 cnt++; 1745 mdiag[i] = 1; 1746 } 1747 } 1748 } 1749 if (!cnt) { 1750 PetscCall(MatShift_Basic(A,v)); 1751 } else { 1752 PetscScalar *olda = a->a; /* preserve pointers to current matrix nonzeros structure and values */ 1753 PetscInt *oldj = a->j, *oldi = a->i; 1754 PetscBool singlemalloc = a->singlemalloc,free_a = a->free_a,free_ij = a->free_ij; 1755 1756 a->a = NULL; 1757 a->j = NULL; 1758 a->i = NULL; 1759 /* increase the values in imax for each row where a diagonal is being inserted then reallocate the matrix data structures */ 1760 for (i=0; i<PetscMin(A->rmap->n,A->cmap->n); i++) { 1761 a->imax[i] += mdiag[i]; 1762 } 1763 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A,0,a->imax)); 1764 1765 /* copy old values into new matrix data structure */ 1766 for (i=0; i<A->rmap->n; i++) { 1767 PetscCall(MatSetValues(A,1,&i,a->imax[i] - mdiag[i],&oldj[oldi[i]],&olda[oldi[i]],ADD_VALUES)); 1768 if (i < A->cmap->n) { 1769 PetscCall(MatSetValue(A,i,i,v,ADD_VALUES)); 1770 } 1771 } 1772 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1773 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1774 if (singlemalloc) { 1775 PetscCall(PetscFree3(olda,oldj,oldi)); 1776 } else { 1777 if (free_a) PetscCall(PetscFree(olda)); 1778 if (free_ij) PetscCall(PetscFree(oldj)); 1779 if (free_ij) PetscCall(PetscFree(oldi)); 1780 } 1781 } 1782 PetscCall(PetscFree(mdiag)); 1783 a->diagonaldense = PETSC_TRUE; 1784 PetscFunctionReturn(0); 1785 } 1786 1787 /* 1788 Checks for missing diagonals 1789 */ 1790 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1791 { 1792 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1793 PetscInt *diag,*ii = a->i,i; 1794 1795 PetscFunctionBegin; 1796 *missing = PETSC_FALSE; 1797 if (A->rmap->n > 0 && !ii) { 1798 *missing = PETSC_TRUE; 1799 if (d) *d = 0; 1800 PetscCall(PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n")); 1801 } else { 1802 PetscInt n; 1803 n = PetscMin(A->rmap->n, A->cmap->n); 1804 diag = a->diag; 1805 for (i=0; i<n; i++) { 1806 if (diag[i] >= ii[i+1]) { 1807 *missing = PETSC_TRUE; 1808 if (d) *d = i; 1809 PetscCall(PetscInfo(A,"Matrix is missing diagonal number %" PetscInt_FMT "\n",i)); 1810 break; 1811 } 1812 } 1813 } 1814 PetscFunctionReturn(0); 1815 } 1816 1817 #include <petscblaslapack.h> 1818 #include <petsc/private/kernels/blockinvert.h> 1819 1820 /* 1821 Note that values is allocated externally by the PC and then passed into this routine 1822 */ 1823 PetscErrorCode MatInvertVariableBlockDiagonal_SeqAIJ(Mat A,PetscInt nblocks,const PetscInt *bsizes,PetscScalar *diag) 1824 { 1825 PetscInt n = A->rmap->n, i, ncnt = 0, *indx,j,bsizemax = 0,*v_pivots; 1826 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 1827 const PetscReal shift = 0.0; 1828 PetscInt ipvt[5]; 1829 PetscScalar work[25],*v_work; 1830 1831 PetscFunctionBegin; 1832 allowzeropivot = PetscNot(A->erroriffailure); 1833 for (i=0; i<nblocks; i++) ncnt += bsizes[i]; 1834 PetscCheck(ncnt == n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Total blocksizes %" PetscInt_FMT " doesn't match number matrix rows %" PetscInt_FMT,ncnt,n); 1835 for (i=0; i<nblocks; i++) { 1836 bsizemax = PetscMax(bsizemax,bsizes[i]); 1837 } 1838 PetscCall(PetscMalloc1(bsizemax,&indx)); 1839 if (bsizemax > 7) { 1840 PetscCall(PetscMalloc2(bsizemax,&v_work,bsizemax,&v_pivots)); 1841 } 1842 ncnt = 0; 1843 for (i=0; i<nblocks; i++) { 1844 for (j=0; j<bsizes[i]; j++) indx[j] = ncnt+j; 1845 PetscCall(MatGetValues(A,bsizes[i],indx,bsizes[i],indx,diag)); 1846 switch (bsizes[i]) { 1847 case 1: 1848 *diag = 1.0/(*diag); 1849 break; 1850 case 2: 1851 PetscCall(PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected)); 1852 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1853 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 1854 break; 1855 case 3: 1856 PetscCall(PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected)); 1857 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1858 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 1859 break; 1860 case 4: 1861 PetscCall(PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected)); 1862 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1863 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 1864 break; 1865 case 5: 1866 PetscCall(PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 1867 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1868 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 1869 break; 1870 case 6: 1871 PetscCall(PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected)); 1872 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1873 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 1874 break; 1875 case 7: 1876 PetscCall(PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected)); 1877 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1878 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 1879 break; 1880 default: 1881 PetscCall(PetscKernel_A_gets_inverse_A(bsizes[i],diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 1882 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1883 PetscCall(PetscKernel_A_gets_transpose_A_N(diag,bsizes[i])); 1884 } 1885 ncnt += bsizes[i]; 1886 diag += bsizes[i]*bsizes[i]; 1887 } 1888 if (bsizemax > 7) { 1889 PetscCall(PetscFree2(v_work,v_pivots)); 1890 } 1891 PetscCall(PetscFree(indx)); 1892 PetscFunctionReturn(0); 1893 } 1894 1895 /* 1896 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1897 */ 1898 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1899 { 1900 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1901 PetscInt i,*diag,m = A->rmap->n; 1902 const MatScalar *v; 1903 PetscScalar *idiag,*mdiag; 1904 1905 PetscFunctionBegin; 1906 if (a->idiagvalid) PetscFunctionReturn(0); 1907 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1908 diag = a->diag; 1909 if (!a->idiag) { 1910 PetscCall(PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work)); 1911 PetscCall(PetscLogObjectMemory((PetscObject)A,3*m*sizeof(PetscScalar))); 1912 } 1913 1914 mdiag = a->mdiag; 1915 idiag = a->idiag; 1916 PetscCall(MatSeqAIJGetArrayRead(A,&v)); 1917 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1918 for (i=0; i<m; i++) { 1919 mdiag[i] = v[diag[i]]; 1920 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1921 if (PetscRealPart(fshift)) { 1922 PetscCall(PetscInfo(A,"Zero diagonal on row %" PetscInt_FMT "\n",i)); 1923 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1924 A->factorerror_zeropivot_value = 0.0; 1925 A->factorerror_zeropivot_row = i; 1926 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %" PetscInt_FMT,i); 1927 } 1928 idiag[i] = 1.0/v[diag[i]]; 1929 } 1930 PetscCall(PetscLogFlops(m)); 1931 } else { 1932 for (i=0; i<m; i++) { 1933 mdiag[i] = v[diag[i]]; 1934 idiag[i] = omega/(fshift + v[diag[i]]); 1935 } 1936 PetscCall(PetscLogFlops(2.0*m)); 1937 } 1938 a->idiagvalid = PETSC_TRUE; 1939 PetscCall(MatSeqAIJRestoreArrayRead(A,&v)); 1940 PetscFunctionReturn(0); 1941 } 1942 1943 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1944 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1945 { 1946 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1947 PetscScalar *x,d,sum,*t,scale; 1948 const MatScalar *v,*idiag=NULL,*mdiag,*aa; 1949 const PetscScalar *b, *bs,*xb, *ts; 1950 PetscInt n,m = A->rmap->n,i; 1951 const PetscInt *idx,*diag; 1952 1953 PetscFunctionBegin; 1954 if (a->inode.use && a->inode.checked && omega == 1.0 && fshift == 0.0) { 1955 PetscCall(MatSOR_SeqAIJ_Inode(A,bb,omega,flag,fshift,its,lits,xx)); 1956 PetscFunctionReturn(0); 1957 } 1958 its = its*lits; 1959 1960 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1961 if (!a->idiagvalid) PetscCall(MatInvertDiagonal_SeqAIJ(A,omega,fshift)); 1962 a->fshift = fshift; 1963 a->omega = omega; 1964 1965 diag = a->diag; 1966 t = a->ssor_work; 1967 idiag = a->idiag; 1968 mdiag = a->mdiag; 1969 1970 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1971 PetscCall(VecGetArray(xx,&x)); 1972 PetscCall(VecGetArrayRead(bb,&b)); 1973 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1974 if (flag == SOR_APPLY_UPPER) { 1975 /* apply (U + D/omega) to the vector */ 1976 bs = b; 1977 for (i=0; i<m; i++) { 1978 d = fshift + mdiag[i]; 1979 n = a->i[i+1] - diag[i] - 1; 1980 idx = a->j + diag[i] + 1; 1981 v = aa + diag[i] + 1; 1982 sum = b[i]*d/omega; 1983 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1984 x[i] = sum; 1985 } 1986 PetscCall(VecRestoreArray(xx,&x)); 1987 PetscCall(VecRestoreArrayRead(bb,&b)); 1988 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1989 PetscCall(PetscLogFlops(a->nz)); 1990 PetscFunctionReturn(0); 1991 } 1992 1993 PetscCheck(flag != SOR_APPLY_LOWER,PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1994 else if (flag & SOR_EISENSTAT) { 1995 /* Let A = L + U + D; where L is lower triangular, 1996 U is upper triangular, E = D/omega; This routine applies 1997 1998 (L + E)^{-1} A (U + E)^{-1} 1999 2000 to a vector efficiently using Eisenstat's trick. 2001 */ 2002 scale = (2.0/omega) - 1.0; 2003 2004 /* x = (E + U)^{-1} b */ 2005 for (i=m-1; i>=0; i--) { 2006 n = a->i[i+1] - diag[i] - 1; 2007 idx = a->j + diag[i] + 1; 2008 v = aa + diag[i] + 1; 2009 sum = b[i]; 2010 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2011 x[i] = sum*idiag[i]; 2012 } 2013 2014 /* t = b - (2*E - D)x */ 2015 v = aa; 2016 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 2017 2018 /* t = (E + L)^{-1}t */ 2019 ts = t; 2020 diag = a->diag; 2021 for (i=0; i<m; i++) { 2022 n = diag[i] - a->i[i]; 2023 idx = a->j + a->i[i]; 2024 v = aa + a->i[i]; 2025 sum = t[i]; 2026 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 2027 t[i] = sum*idiag[i]; 2028 /* x = x + t */ 2029 x[i] += t[i]; 2030 } 2031 2032 PetscCall(PetscLogFlops(6.0*m-1 + 2.0*a->nz)); 2033 PetscCall(VecRestoreArray(xx,&x)); 2034 PetscCall(VecRestoreArrayRead(bb,&b)); 2035 PetscFunctionReturn(0); 2036 } 2037 if (flag & SOR_ZERO_INITIAL_GUESS) { 2038 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2039 for (i=0; i<m; i++) { 2040 n = diag[i] - a->i[i]; 2041 idx = a->j + a->i[i]; 2042 v = aa + a->i[i]; 2043 sum = b[i]; 2044 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2045 t[i] = sum; 2046 x[i] = sum*idiag[i]; 2047 } 2048 xb = t; 2049 PetscCall(PetscLogFlops(a->nz)); 2050 } else xb = b; 2051 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2052 for (i=m-1; i>=0; i--) { 2053 n = a->i[i+1] - diag[i] - 1; 2054 idx = a->j + diag[i] + 1; 2055 v = aa + diag[i] + 1; 2056 sum = xb[i]; 2057 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2058 if (xb == b) { 2059 x[i] = sum*idiag[i]; 2060 } else { 2061 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2062 } 2063 } 2064 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2065 } 2066 its--; 2067 } 2068 while (its--) { 2069 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2070 for (i=0; i<m; i++) { 2071 /* lower */ 2072 n = diag[i] - a->i[i]; 2073 idx = a->j + a->i[i]; 2074 v = aa + a->i[i]; 2075 sum = b[i]; 2076 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2077 t[i] = sum; /* save application of the lower-triangular part */ 2078 /* upper */ 2079 n = a->i[i+1] - diag[i] - 1; 2080 idx = a->j + diag[i] + 1; 2081 v = aa + diag[i] + 1; 2082 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2083 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2084 } 2085 xb = t; 2086 PetscCall(PetscLogFlops(2.0*a->nz)); 2087 } else xb = b; 2088 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2089 for (i=m-1; i>=0; i--) { 2090 sum = xb[i]; 2091 if (xb == b) { 2092 /* whole matrix (no checkpointing available) */ 2093 n = a->i[i+1] - a->i[i]; 2094 idx = a->j + a->i[i]; 2095 v = aa + a->i[i]; 2096 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2097 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 2098 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 2099 n = a->i[i+1] - diag[i] - 1; 2100 idx = a->j + diag[i] + 1; 2101 v = aa + diag[i] + 1; 2102 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2103 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2104 } 2105 } 2106 if (xb == b) { 2107 PetscCall(PetscLogFlops(2.0*a->nz)); 2108 } else { 2109 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2110 } 2111 } 2112 } 2113 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2114 PetscCall(VecRestoreArray(xx,&x)); 2115 PetscCall(VecRestoreArrayRead(bb,&b)); 2116 PetscFunctionReturn(0); 2117 } 2118 2119 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 2120 { 2121 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2122 2123 PetscFunctionBegin; 2124 info->block_size = 1.0; 2125 info->nz_allocated = a->maxnz; 2126 info->nz_used = a->nz; 2127 info->nz_unneeded = (a->maxnz - a->nz); 2128 info->assemblies = A->num_ass; 2129 info->mallocs = A->info.mallocs; 2130 info->memory = ((PetscObject)A)->mem; 2131 if (A->factortype) { 2132 info->fill_ratio_given = A->info.fill_ratio_given; 2133 info->fill_ratio_needed = A->info.fill_ratio_needed; 2134 info->factor_mallocs = A->info.factor_mallocs; 2135 } else { 2136 info->fill_ratio_given = 0; 2137 info->fill_ratio_needed = 0; 2138 info->factor_mallocs = 0; 2139 } 2140 PetscFunctionReturn(0); 2141 } 2142 2143 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2144 { 2145 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2146 PetscInt i,m = A->rmap->n - 1; 2147 const PetscScalar *xx; 2148 PetscScalar *bb,*aa; 2149 PetscInt d = 0; 2150 2151 PetscFunctionBegin; 2152 if (x && b) { 2153 PetscCall(VecGetArrayRead(x,&xx)); 2154 PetscCall(VecGetArray(b,&bb)); 2155 for (i=0; i<N; i++) { 2156 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2157 if (rows[i] >= A->cmap->n) continue; 2158 bb[rows[i]] = diag*xx[rows[i]]; 2159 } 2160 PetscCall(VecRestoreArrayRead(x,&xx)); 2161 PetscCall(VecRestoreArray(b,&bb)); 2162 } 2163 2164 PetscCall(MatSeqAIJGetArray(A,&aa)); 2165 if (a->keepnonzeropattern) { 2166 for (i=0; i<N; i++) { 2167 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2168 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]])); 2169 } 2170 if (diag != 0.0) { 2171 for (i=0; i<N; i++) { 2172 d = rows[i]; 2173 if (rows[i] >= A->cmap->n) continue; 2174 PetscCheck(a->diag[d] < a->i[d+1],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in the zeroed row %" PetscInt_FMT,d); 2175 } 2176 for (i=0; i<N; i++) { 2177 if (rows[i] >= A->cmap->n) continue; 2178 aa[a->diag[rows[i]]] = diag; 2179 } 2180 } 2181 } else { 2182 if (diag != 0.0) { 2183 for (i=0; i<N; i++) { 2184 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2185 if (a->ilen[rows[i]] > 0) { 2186 if (rows[i] >= A->cmap->n) { 2187 a->ilen[rows[i]] = 0; 2188 } else { 2189 a->ilen[rows[i]] = 1; 2190 aa[a->i[rows[i]]] = diag; 2191 a->j[a->i[rows[i]]] = rows[i]; 2192 } 2193 } else if (rows[i] < A->cmap->n) { /* in case row was completely empty */ 2194 PetscCall(MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES)); 2195 } 2196 } 2197 } else { 2198 for (i=0; i<N; i++) { 2199 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2200 a->ilen[rows[i]] = 0; 2201 } 2202 } 2203 A->nonzerostate++; 2204 } 2205 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 2206 PetscCall((*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY)); 2207 PetscFunctionReturn(0); 2208 } 2209 2210 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2211 { 2212 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2213 PetscInt i,j,m = A->rmap->n - 1,d = 0; 2214 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 2215 const PetscScalar *xx; 2216 PetscScalar *bb,*aa; 2217 2218 PetscFunctionBegin; 2219 if (!N) PetscFunctionReturn(0); 2220 PetscCall(MatSeqAIJGetArray(A,&aa)); 2221 if (x && b) { 2222 PetscCall(VecGetArrayRead(x,&xx)); 2223 PetscCall(VecGetArray(b,&bb)); 2224 vecs = PETSC_TRUE; 2225 } 2226 PetscCall(PetscCalloc1(A->rmap->n,&zeroed)); 2227 for (i=0; i<N; i++) { 2228 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2229 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]])); 2230 2231 zeroed[rows[i]] = PETSC_TRUE; 2232 } 2233 for (i=0; i<A->rmap->n; i++) { 2234 if (!zeroed[i]) { 2235 for (j=a->i[i]; j<a->i[i+1]; j++) { 2236 if (a->j[j] < A->rmap->n && zeroed[a->j[j]]) { 2237 if (vecs) bb[i] -= aa[j]*xx[a->j[j]]; 2238 aa[j] = 0.0; 2239 } 2240 } 2241 } else if (vecs && i < A->cmap->N) bb[i] = diag*xx[i]; 2242 } 2243 if (x && b) { 2244 PetscCall(VecRestoreArrayRead(x,&xx)); 2245 PetscCall(VecRestoreArray(b,&bb)); 2246 } 2247 PetscCall(PetscFree(zeroed)); 2248 if (diag != 0.0) { 2249 PetscCall(MatMissingDiagonal_SeqAIJ(A,&missing,&d)); 2250 if (missing) { 2251 for (i=0; i<N; i++) { 2252 if (rows[i] >= A->cmap->N) continue; 2253 PetscCheck(!a->nonew || rows[i] < d,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %" PetscInt_FMT " (%" PetscInt_FMT ")",d,rows[i]); 2254 PetscCall(MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES)); 2255 } 2256 } else { 2257 for (i=0; i<N; i++) { 2258 aa[a->diag[rows[i]]] = diag; 2259 } 2260 } 2261 } 2262 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 2263 PetscCall((*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY)); 2264 PetscFunctionReturn(0); 2265 } 2266 2267 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2268 { 2269 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2270 const PetscScalar *aa; 2271 PetscInt *itmp; 2272 2273 PetscFunctionBegin; 2274 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2275 *nz = a->i[row+1] - a->i[row]; 2276 if (v) *v = (PetscScalar*)(aa + a->i[row]); 2277 if (idx) { 2278 itmp = a->j + a->i[row]; 2279 if (*nz) *idx = itmp; 2280 else *idx = NULL; 2281 } 2282 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2283 PetscFunctionReturn(0); 2284 } 2285 2286 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2287 { 2288 PetscFunctionBegin; 2289 if (nz) *nz = 0; 2290 if (idx) *idx = NULL; 2291 if (v) *v = NULL; 2292 PetscFunctionReturn(0); 2293 } 2294 2295 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 2296 { 2297 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2298 const MatScalar *v; 2299 PetscReal sum = 0.0; 2300 PetscInt i,j; 2301 2302 PetscFunctionBegin; 2303 PetscCall(MatSeqAIJGetArrayRead(A,&v)); 2304 if (type == NORM_FROBENIUS) { 2305 #if defined(PETSC_USE_REAL___FP16) 2306 PetscBLASInt one = 1,nz = a->nz; 2307 PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&nz,v,&one)); 2308 #else 2309 for (i=0; i<a->nz; i++) { 2310 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2311 } 2312 *nrm = PetscSqrtReal(sum); 2313 #endif 2314 PetscCall(PetscLogFlops(2.0*a->nz)); 2315 } else if (type == NORM_1) { 2316 PetscReal *tmp; 2317 PetscInt *jj = a->j; 2318 PetscCall(PetscCalloc1(A->cmap->n+1,&tmp)); 2319 *nrm = 0.0; 2320 for (j=0; j<a->nz; j++) { 2321 tmp[*jj++] += PetscAbsScalar(*v); v++; 2322 } 2323 for (j=0; j<A->cmap->n; j++) { 2324 if (tmp[j] > *nrm) *nrm = tmp[j]; 2325 } 2326 PetscCall(PetscFree(tmp)); 2327 PetscCall(PetscLogFlops(PetscMax(a->nz-1,0))); 2328 } else if (type == NORM_INFINITY) { 2329 *nrm = 0.0; 2330 for (j=0; j<A->rmap->n; j++) { 2331 const PetscScalar *v2 = v + a->i[j]; 2332 sum = 0.0; 2333 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2334 sum += PetscAbsScalar(*v2); v2++; 2335 } 2336 if (sum > *nrm) *nrm = sum; 2337 } 2338 PetscCall(PetscLogFlops(PetscMax(a->nz-1,0))); 2339 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2340 PetscCall(MatSeqAIJRestoreArrayRead(A,&v)); 2341 PetscFunctionReturn(0); 2342 } 2343 2344 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 2345 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 2346 { 2347 PetscInt i,j,anzj; 2348 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 2349 PetscInt an=A->cmap->N,am=A->rmap->N; 2350 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2351 2352 PetscFunctionBegin; 2353 /* Allocate space for symbolic transpose info and work array */ 2354 PetscCall(PetscCalloc1(an+1,&ati)); 2355 PetscCall(PetscMalloc1(ai[am],&atj)); 2356 PetscCall(PetscMalloc1(an,&atfill)); 2357 2358 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2359 /* Note: offset by 1 for fast conversion into csr format. */ 2360 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2361 /* Form ati for csr format of A^T. */ 2362 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2363 2364 /* Copy ati into atfill so we have locations of the next free space in atj */ 2365 PetscCall(PetscArraycpy(atfill,ati,an)); 2366 2367 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2368 for (i=0;i<am;i++) { 2369 anzj = ai[i+1] - ai[i]; 2370 for (j=0;j<anzj;j++) { 2371 atj[atfill[*aj]] = i; 2372 atfill[*aj++] += 1; 2373 } 2374 } 2375 2376 /* Clean up temporary space and complete requests. */ 2377 PetscCall(PetscFree(atfill)); 2378 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B)); 2379 PetscCall(MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs))); 2380 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 2381 2382 b = (Mat_SeqAIJ*)((*B)->data); 2383 b->free_a = PETSC_FALSE; 2384 b->free_ij = PETSC_TRUE; 2385 b->nonew = 0; 2386 PetscFunctionReturn(0); 2387 } 2388 2389 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2390 { 2391 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2392 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2393 const MatScalar *va,*vb; 2394 PetscInt ma,na,mb,nb, i; 2395 2396 PetscFunctionBegin; 2397 PetscCall(MatGetSize(A,&ma,&na)); 2398 PetscCall(MatGetSize(B,&mb,&nb)); 2399 if (ma!=nb || na!=mb) { 2400 *f = PETSC_FALSE; 2401 PetscFunctionReturn(0); 2402 } 2403 PetscCall(MatSeqAIJGetArrayRead(A,&va)); 2404 PetscCall(MatSeqAIJGetArrayRead(B,&vb)); 2405 aii = aij->i; bii = bij->i; 2406 adx = aij->j; bdx = bij->j; 2407 PetscCall(PetscMalloc1(ma,&aptr)); 2408 PetscCall(PetscMalloc1(mb,&bptr)); 2409 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2410 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2411 2412 *f = PETSC_TRUE; 2413 for (i=0; i<ma; i++) { 2414 while (aptr[i]<aii[i+1]) { 2415 PetscInt idc,idr; 2416 PetscScalar vc,vr; 2417 /* column/row index/value */ 2418 idc = adx[aptr[i]]; 2419 idr = bdx[bptr[idc]]; 2420 vc = va[aptr[i]]; 2421 vr = vb[bptr[idc]]; 2422 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2423 *f = PETSC_FALSE; 2424 goto done; 2425 } else { 2426 aptr[i]++; 2427 if (B || i!=idc) bptr[idc]++; 2428 } 2429 } 2430 } 2431 done: 2432 PetscCall(PetscFree(aptr)); 2433 PetscCall(PetscFree(bptr)); 2434 PetscCall(MatSeqAIJRestoreArrayRead(A,&va)); 2435 PetscCall(MatSeqAIJRestoreArrayRead(B,&vb)); 2436 PetscFunctionReturn(0); 2437 } 2438 2439 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2440 { 2441 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2442 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2443 MatScalar *va,*vb; 2444 PetscInt ma,na,mb,nb, i; 2445 2446 PetscFunctionBegin; 2447 PetscCall(MatGetSize(A,&ma,&na)); 2448 PetscCall(MatGetSize(B,&mb,&nb)); 2449 if (ma!=nb || na!=mb) { 2450 *f = PETSC_FALSE; 2451 PetscFunctionReturn(0); 2452 } 2453 aii = aij->i; bii = bij->i; 2454 adx = aij->j; bdx = bij->j; 2455 va = aij->a; vb = bij->a; 2456 PetscCall(PetscMalloc1(ma,&aptr)); 2457 PetscCall(PetscMalloc1(mb,&bptr)); 2458 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2459 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2460 2461 *f = PETSC_TRUE; 2462 for (i=0; i<ma; i++) { 2463 while (aptr[i]<aii[i+1]) { 2464 PetscInt idc,idr; 2465 PetscScalar vc,vr; 2466 /* column/row index/value */ 2467 idc = adx[aptr[i]]; 2468 idr = bdx[bptr[idc]]; 2469 vc = va[aptr[i]]; 2470 vr = vb[bptr[idc]]; 2471 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2472 *f = PETSC_FALSE; 2473 goto done; 2474 } else { 2475 aptr[i]++; 2476 if (B || i!=idc) bptr[idc]++; 2477 } 2478 } 2479 } 2480 done: 2481 PetscCall(PetscFree(aptr)); 2482 PetscCall(PetscFree(bptr)); 2483 PetscFunctionReturn(0); 2484 } 2485 2486 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2487 { 2488 PetscFunctionBegin; 2489 PetscCall(MatIsTranspose_SeqAIJ(A,A,tol,f)); 2490 PetscFunctionReturn(0); 2491 } 2492 2493 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2494 { 2495 PetscFunctionBegin; 2496 PetscCall(MatIsHermitianTranspose_SeqAIJ(A,A,tol,f)); 2497 PetscFunctionReturn(0); 2498 } 2499 2500 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2501 { 2502 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2503 const PetscScalar *l,*r; 2504 PetscScalar x; 2505 MatScalar *v; 2506 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz; 2507 const PetscInt *jj; 2508 2509 PetscFunctionBegin; 2510 if (ll) { 2511 /* The local size is used so that VecMPI can be passed to this routine 2512 by MatDiagonalScale_MPIAIJ */ 2513 PetscCall(VecGetLocalSize(ll,&m)); 2514 PetscCheck(m == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2515 PetscCall(VecGetArrayRead(ll,&l)); 2516 PetscCall(MatSeqAIJGetArray(A,&v)); 2517 for (i=0; i<m; i++) { 2518 x = l[i]; 2519 M = a->i[i+1] - a->i[i]; 2520 for (j=0; j<M; j++) (*v++) *= x; 2521 } 2522 PetscCall(VecRestoreArrayRead(ll,&l)); 2523 PetscCall(PetscLogFlops(nz)); 2524 PetscCall(MatSeqAIJRestoreArray(A,&v)); 2525 } 2526 if (rr) { 2527 PetscCall(VecGetLocalSize(rr,&n)); 2528 PetscCheck(n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2529 PetscCall(VecGetArrayRead(rr,&r)); 2530 PetscCall(MatSeqAIJGetArray(A,&v)); 2531 jj = a->j; 2532 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2533 PetscCall(MatSeqAIJRestoreArray(A,&v)); 2534 PetscCall(VecRestoreArrayRead(rr,&r)); 2535 PetscCall(PetscLogFlops(nz)); 2536 } 2537 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 2538 PetscFunctionReturn(0); 2539 } 2540 2541 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2542 { 2543 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2544 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2545 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2546 const PetscInt *irow,*icol; 2547 const PetscScalar *aa; 2548 PetscInt nrows,ncols; 2549 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2550 MatScalar *a_new,*mat_a; 2551 Mat C; 2552 PetscBool stride; 2553 2554 PetscFunctionBegin; 2555 PetscCall(ISGetIndices(isrow,&irow)); 2556 PetscCall(ISGetLocalSize(isrow,&nrows)); 2557 PetscCall(ISGetLocalSize(iscol,&ncols)); 2558 2559 PetscCall(PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride)); 2560 if (stride) { 2561 PetscCall(ISStrideGetInfo(iscol,&first,&step)); 2562 } else { 2563 first = 0; 2564 step = 0; 2565 } 2566 if (stride && step == 1) { 2567 /* special case of contiguous rows */ 2568 PetscCall(PetscMalloc2(nrows,&lens,nrows,&starts)); 2569 /* loop over new rows determining lens and starting points */ 2570 for (i=0; i<nrows; i++) { 2571 kstart = ai[irow[i]]; 2572 kend = kstart + ailen[irow[i]]; 2573 starts[i] = kstart; 2574 for (k=kstart; k<kend; k++) { 2575 if (aj[k] >= first) { 2576 starts[i] = k; 2577 break; 2578 } 2579 } 2580 sum = 0; 2581 while (k < kend) { 2582 if (aj[k++] >= first+ncols) break; 2583 sum++; 2584 } 2585 lens[i] = sum; 2586 } 2587 /* create submatrix */ 2588 if (scall == MAT_REUSE_MATRIX) { 2589 PetscInt n_cols,n_rows; 2590 PetscCall(MatGetSize(*B,&n_rows,&n_cols)); 2591 PetscCheck(n_rows == nrows && n_cols == ncols,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2592 PetscCall(MatZeroEntries(*B)); 2593 C = *B; 2594 } else { 2595 PetscInt rbs,cbs; 2596 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2597 PetscCall(MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE)); 2598 PetscCall(ISGetBlockSize(isrow,&rbs)); 2599 PetscCall(ISGetBlockSize(iscol,&cbs)); 2600 PetscCall(MatSetBlockSizes(C,rbs,cbs)); 2601 PetscCall(MatSetType(C,((PetscObject)A)->type_name)); 2602 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens)); 2603 } 2604 c = (Mat_SeqAIJ*)C->data; 2605 2606 /* loop over rows inserting into submatrix */ 2607 a_new = c->a; 2608 j_new = c->j; 2609 i_new = c->i; 2610 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2611 for (i=0; i<nrows; i++) { 2612 ii = starts[i]; 2613 lensi = lens[i]; 2614 for (k=0; k<lensi; k++) { 2615 *j_new++ = aj[ii+k] - first; 2616 } 2617 PetscCall(PetscArraycpy(a_new,aa + starts[i],lensi)); 2618 a_new += lensi; 2619 i_new[i+1] = i_new[i] + lensi; 2620 c->ilen[i] = lensi; 2621 } 2622 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2623 PetscCall(PetscFree2(lens,starts)); 2624 } else { 2625 PetscCall(ISGetIndices(iscol,&icol)); 2626 PetscCall(PetscCalloc1(oldcols,&smap)); 2627 PetscCall(PetscMalloc1(1+nrows,&lens)); 2628 for (i=0; i<ncols; i++) { 2629 PetscCheck(icol[i] < oldcols,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requesting column beyond largest column icol[%" PetscInt_FMT "] %" PetscInt_FMT " >= A->cmap->n %" PetscInt_FMT,i,icol[i],oldcols); 2630 smap[icol[i]] = i+1; 2631 } 2632 2633 /* determine lens of each row */ 2634 for (i=0; i<nrows; i++) { 2635 kstart = ai[irow[i]]; 2636 kend = kstart + a->ilen[irow[i]]; 2637 lens[i] = 0; 2638 for (k=kstart; k<kend; k++) { 2639 if (smap[aj[k]]) { 2640 lens[i]++; 2641 } 2642 } 2643 } 2644 /* Create and fill new matrix */ 2645 if (scall == MAT_REUSE_MATRIX) { 2646 PetscBool equal; 2647 2648 c = (Mat_SeqAIJ*)((*B)->data); 2649 PetscCheck((*B)->rmap->n == nrows && (*B)->cmap->n == ncols,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2650 PetscCall(PetscArraycmp(c->ilen,lens,(*B)->rmap->n,&equal)); 2651 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2652 PetscCall(PetscArrayzero(c->ilen,(*B)->rmap->n)); 2653 C = *B; 2654 } else { 2655 PetscInt rbs,cbs; 2656 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2657 PetscCall(MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE)); 2658 PetscCall(ISGetBlockSize(isrow,&rbs)); 2659 PetscCall(ISGetBlockSize(iscol,&cbs)); 2660 PetscCall(MatSetBlockSizes(C,rbs,cbs)); 2661 PetscCall(MatSetType(C,((PetscObject)A)->type_name)); 2662 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens)); 2663 } 2664 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2665 c = (Mat_SeqAIJ*)(C->data); 2666 for (i=0; i<nrows; i++) { 2667 row = irow[i]; 2668 kstart = ai[row]; 2669 kend = kstart + a->ilen[row]; 2670 mat_i = c->i[i]; 2671 mat_j = c->j + mat_i; 2672 mat_a = c->a + mat_i; 2673 mat_ilen = c->ilen + i; 2674 for (k=kstart; k<kend; k++) { 2675 if ((tcol=smap[a->j[k]])) { 2676 *mat_j++ = tcol - 1; 2677 *mat_a++ = aa[k]; 2678 (*mat_ilen)++; 2679 2680 } 2681 } 2682 } 2683 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2684 /* Free work space */ 2685 PetscCall(ISRestoreIndices(iscol,&icol)); 2686 PetscCall(PetscFree(smap)); 2687 PetscCall(PetscFree(lens)); 2688 /* sort */ 2689 for (i = 0; i < nrows; i++) { 2690 PetscInt ilen; 2691 2692 mat_i = c->i[i]; 2693 mat_j = c->j + mat_i; 2694 mat_a = c->a + mat_i; 2695 ilen = c->ilen[i]; 2696 PetscCall(PetscSortIntWithScalarArray(ilen,mat_j,mat_a)); 2697 } 2698 } 2699 #if defined(PETSC_HAVE_DEVICE) 2700 PetscCall(MatBindToCPU(C,A->boundtocpu)); 2701 #endif 2702 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 2703 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 2704 2705 PetscCall(ISRestoreIndices(isrow,&irow)); 2706 *B = C; 2707 PetscFunctionReturn(0); 2708 } 2709 2710 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2711 { 2712 Mat B; 2713 2714 PetscFunctionBegin; 2715 if (scall == MAT_INITIAL_MATRIX) { 2716 PetscCall(MatCreate(subComm,&B)); 2717 PetscCall(MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n)); 2718 PetscCall(MatSetBlockSizesFromMats(B,mat,mat)); 2719 PetscCall(MatSetType(B,MATSEQAIJ)); 2720 PetscCall(MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE)); 2721 *subMat = B; 2722 } else { 2723 PetscCall(MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN)); 2724 } 2725 PetscFunctionReturn(0); 2726 } 2727 2728 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2729 { 2730 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2731 Mat outA; 2732 PetscBool row_identity,col_identity; 2733 2734 PetscFunctionBegin; 2735 PetscCheck(info->levels == 0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2736 2737 PetscCall(ISIdentity(row,&row_identity)); 2738 PetscCall(ISIdentity(col,&col_identity)); 2739 2740 outA = inA; 2741 outA->factortype = MAT_FACTOR_LU; 2742 PetscCall(PetscFree(inA->solvertype)); 2743 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype)); 2744 2745 PetscCall(PetscObjectReference((PetscObject)row)); 2746 PetscCall(ISDestroy(&a->row)); 2747 2748 a->row = row; 2749 2750 PetscCall(PetscObjectReference((PetscObject)col)); 2751 PetscCall(ISDestroy(&a->col)); 2752 2753 a->col = col; 2754 2755 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2756 PetscCall(ISDestroy(&a->icol)); 2757 PetscCall(ISInvertPermutation(col,PETSC_DECIDE,&a->icol)); 2758 PetscCall(PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol)); 2759 2760 if (!a->solve_work) { /* this matrix may have been factored before */ 2761 PetscCall(PetscMalloc1(inA->rmap->n+1,&a->solve_work)); 2762 PetscCall(PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar))); 2763 } 2764 2765 PetscCall(MatMarkDiagonal_SeqAIJ(inA)); 2766 if (row_identity && col_identity) { 2767 PetscCall(MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info)); 2768 } else { 2769 PetscCall(MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info)); 2770 } 2771 PetscFunctionReturn(0); 2772 } 2773 2774 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2775 { 2776 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2777 PetscScalar *v; 2778 PetscBLASInt one = 1,bnz; 2779 2780 PetscFunctionBegin; 2781 PetscCall(MatSeqAIJGetArray(inA,&v)); 2782 PetscCall(PetscBLASIntCast(a->nz,&bnz)); 2783 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&alpha,v,&one)); 2784 PetscCall(PetscLogFlops(a->nz)); 2785 PetscCall(MatSeqAIJRestoreArray(inA,&v)); 2786 PetscCall(MatSeqAIJInvalidateDiagonal(inA)); 2787 PetscFunctionReturn(0); 2788 } 2789 2790 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2791 { 2792 PetscInt i; 2793 2794 PetscFunctionBegin; 2795 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2796 PetscCall(PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr)); 2797 2798 for (i=0; i<submatj->nrqr; ++i) { 2799 PetscCall(PetscFree(submatj->sbuf2[i])); 2800 } 2801 PetscCall(PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1)); 2802 2803 if (submatj->rbuf1) { 2804 PetscCall(PetscFree(submatj->rbuf1[0])); 2805 PetscCall(PetscFree(submatj->rbuf1)); 2806 } 2807 2808 for (i=0; i<submatj->nrqs; ++i) { 2809 PetscCall(PetscFree(submatj->rbuf3[i])); 2810 } 2811 PetscCall(PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3)); 2812 PetscCall(PetscFree(submatj->pa)); 2813 } 2814 2815 #if defined(PETSC_USE_CTABLE) 2816 PetscCall(PetscTableDestroy((PetscTable*)&submatj->rmap)); 2817 if (submatj->cmap_loc) PetscCall(PetscFree(submatj->cmap_loc)); 2818 PetscCall(PetscFree(submatj->rmap_loc)); 2819 #else 2820 PetscCall(PetscFree(submatj->rmap)); 2821 #endif 2822 2823 if (!submatj->allcolumns) { 2824 #if defined(PETSC_USE_CTABLE) 2825 PetscCall(PetscTableDestroy((PetscTable*)&submatj->cmap)); 2826 #else 2827 PetscCall(PetscFree(submatj->cmap)); 2828 #endif 2829 } 2830 PetscCall(PetscFree(submatj->row2proc)); 2831 2832 PetscCall(PetscFree(submatj)); 2833 PetscFunctionReturn(0); 2834 } 2835 2836 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2837 { 2838 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2839 Mat_SubSppt *submatj = c->submatis1; 2840 2841 PetscFunctionBegin; 2842 PetscCall((*submatj->destroy)(C)); 2843 PetscCall(MatDestroySubMatrix_Private(submatj)); 2844 PetscFunctionReturn(0); 2845 } 2846 2847 /* Note this has code duplication with MatDestroySubMatrices_SeqBAIJ() */ 2848 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[]) 2849 { 2850 PetscInt i; 2851 Mat C; 2852 Mat_SeqAIJ *c; 2853 Mat_SubSppt *submatj; 2854 2855 PetscFunctionBegin; 2856 for (i=0; i<n; i++) { 2857 C = (*mat)[i]; 2858 c = (Mat_SeqAIJ*)C->data; 2859 submatj = c->submatis1; 2860 if (submatj) { 2861 if (--((PetscObject)C)->refct <= 0) { 2862 PetscCall(PetscFree(C->factorprefix)); 2863 PetscCall((*submatj->destroy)(C)); 2864 PetscCall(MatDestroySubMatrix_Private(submatj)); 2865 PetscCall(PetscFree(C->defaultvectype)); 2866 PetscCall(PetscLayoutDestroy(&C->rmap)); 2867 PetscCall(PetscLayoutDestroy(&C->cmap)); 2868 PetscCall(PetscHeaderDestroy(&C)); 2869 } 2870 } else { 2871 PetscCall(MatDestroy(&C)); 2872 } 2873 } 2874 2875 /* Destroy Dummy submatrices created for reuse */ 2876 PetscCall(MatDestroySubMatrices_Dummy(n,mat)); 2877 2878 PetscCall(PetscFree(*mat)); 2879 PetscFunctionReturn(0); 2880 } 2881 2882 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2883 { 2884 PetscInt i; 2885 2886 PetscFunctionBegin; 2887 if (scall == MAT_INITIAL_MATRIX) { 2888 PetscCall(PetscCalloc1(n+1,B)); 2889 } 2890 2891 for (i=0; i<n; i++) { 2892 PetscCall(MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i])); 2893 } 2894 PetscFunctionReturn(0); 2895 } 2896 2897 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2898 { 2899 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2900 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2901 const PetscInt *idx; 2902 PetscInt start,end,*ai,*aj; 2903 PetscBT table; 2904 2905 PetscFunctionBegin; 2906 m = A->rmap->n; 2907 ai = a->i; 2908 aj = a->j; 2909 2910 PetscCheck(ov >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2911 2912 PetscCall(PetscMalloc1(m+1,&nidx)); 2913 PetscCall(PetscBTCreate(m,&table)); 2914 2915 for (i=0; i<is_max; i++) { 2916 /* Initialize the two local arrays */ 2917 isz = 0; 2918 PetscCall(PetscBTMemzero(m,table)); 2919 2920 /* Extract the indices, assume there can be duplicate entries */ 2921 PetscCall(ISGetIndices(is[i],&idx)); 2922 PetscCall(ISGetLocalSize(is[i],&n)); 2923 2924 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2925 for (j=0; j<n; ++j) { 2926 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2927 } 2928 PetscCall(ISRestoreIndices(is[i],&idx)); 2929 PetscCall(ISDestroy(&is[i])); 2930 2931 k = 0; 2932 for (j=0; j<ov; j++) { /* for each overlap */ 2933 n = isz; 2934 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2935 row = nidx[k]; 2936 start = ai[row]; 2937 end = ai[row+1]; 2938 for (l = start; l<end; l++) { 2939 val = aj[l]; 2940 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2941 } 2942 } 2943 } 2944 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i))); 2945 } 2946 PetscCall(PetscBTDestroy(&table)); 2947 PetscCall(PetscFree(nidx)); 2948 PetscFunctionReturn(0); 2949 } 2950 2951 /* -------------------------------------------------------------- */ 2952 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2953 { 2954 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2955 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2956 const PetscInt *row,*col; 2957 PetscInt *cnew,j,*lens; 2958 IS icolp,irowp; 2959 PetscInt *cwork = NULL; 2960 PetscScalar *vwork = NULL; 2961 2962 PetscFunctionBegin; 2963 PetscCall(ISInvertPermutation(rowp,PETSC_DECIDE,&irowp)); 2964 PetscCall(ISGetIndices(irowp,&row)); 2965 PetscCall(ISInvertPermutation(colp,PETSC_DECIDE,&icolp)); 2966 PetscCall(ISGetIndices(icolp,&col)); 2967 2968 /* determine lengths of permuted rows */ 2969 PetscCall(PetscMalloc1(m+1,&lens)); 2970 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2971 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 2972 PetscCall(MatSetSizes(*B,m,n,m,n)); 2973 PetscCall(MatSetBlockSizesFromMats(*B,A,A)); 2974 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 2975 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens)); 2976 PetscCall(PetscFree(lens)); 2977 2978 PetscCall(PetscMalloc1(n,&cnew)); 2979 for (i=0; i<m; i++) { 2980 PetscCall(MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork)); 2981 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2982 PetscCall(MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES)); 2983 PetscCall(MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork)); 2984 } 2985 PetscCall(PetscFree(cnew)); 2986 2987 (*B)->assembled = PETSC_FALSE; 2988 2989 #if defined(PETSC_HAVE_DEVICE) 2990 PetscCall(MatBindToCPU(*B,A->boundtocpu)); 2991 #endif 2992 PetscCall(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY)); 2993 PetscCall(MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY)); 2994 PetscCall(ISRestoreIndices(irowp,&row)); 2995 PetscCall(ISRestoreIndices(icolp,&col)); 2996 PetscCall(ISDestroy(&irowp)); 2997 PetscCall(ISDestroy(&icolp)); 2998 if (rowp == colp) { 2999 PetscCall(MatPropagateSymmetryOptions(A,*B)); 3000 } 3001 PetscFunctionReturn(0); 3002 } 3003 3004 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 3005 { 3006 PetscFunctionBegin; 3007 /* If the two matrices have the same copy implementation, use fast copy. */ 3008 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 3009 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3010 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 3011 const PetscScalar *aa; 3012 3013 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 3014 PetscCheck(a->i[A->rmap->n] == b->i[B->rmap->n],PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different %" PetscInt_FMT " != %" PetscInt_FMT,a->i[A->rmap->n],b->i[B->rmap->n]); 3015 PetscCall(PetscArraycpy(b->a,aa,a->i[A->rmap->n])); 3016 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 3017 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 3018 } else { 3019 PetscCall(MatCopy_Basic(A,B,str)); 3020 } 3021 PetscFunctionReturn(0); 3022 } 3023 3024 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 3025 { 3026 PetscFunctionBegin; 3027 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,NULL)); 3028 PetscFunctionReturn(0); 3029 } 3030 3031 PETSC_INTERN PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 3032 { 3033 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3034 3035 PetscFunctionBegin; 3036 *array = a->a; 3037 PetscFunctionReturn(0); 3038 } 3039 3040 PETSC_INTERN PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 3041 { 3042 PetscFunctionBegin; 3043 *array = NULL; 3044 PetscFunctionReturn(0); 3045 } 3046 3047 /* 3048 Computes the number of nonzeros per row needed for preallocation when X and Y 3049 have different nonzero structure. 3050 */ 3051 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 3052 { 3053 PetscInt i,j,k,nzx,nzy; 3054 3055 PetscFunctionBegin; 3056 /* Set the number of nonzeros in the new matrix */ 3057 for (i=0; i<m; i++) { 3058 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 3059 nzx = xi[i+1] - xi[i]; 3060 nzy = yi[i+1] - yi[i]; 3061 nnz[i] = 0; 3062 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 3063 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 3064 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 3065 nnz[i]++; 3066 } 3067 for (; k<nzy; k++) nnz[i]++; 3068 } 3069 PetscFunctionReturn(0); 3070 } 3071 3072 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 3073 { 3074 PetscInt m = Y->rmap->N; 3075 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 3076 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 3077 3078 PetscFunctionBegin; 3079 /* Set the number of nonzeros in the new matrix */ 3080 PetscCall(MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz)); 3081 PetscFunctionReturn(0); 3082 } 3083 3084 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 3085 { 3086 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 3087 3088 PetscFunctionBegin; 3089 if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) { 3090 PetscBool e = x->nz == y->nz ? PETSC_TRUE : PETSC_FALSE; 3091 if (e) { 3092 PetscCall(PetscArraycmp(x->i,y->i,Y->rmap->n+1,&e)); 3093 if (e) { 3094 PetscCall(PetscArraycmp(x->j,y->j,y->nz,&e)); 3095 if (e) str = SAME_NONZERO_PATTERN; 3096 } 3097 } 3098 if (!e) PetscCheck(str != SAME_NONZERO_PATTERN,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatStructure is not SAME_NONZERO_PATTERN"); 3099 } 3100 if (str == SAME_NONZERO_PATTERN) { 3101 const PetscScalar *xa; 3102 PetscScalar *ya,alpha = a; 3103 PetscBLASInt one = 1,bnz; 3104 3105 PetscCall(PetscBLASIntCast(x->nz,&bnz)); 3106 PetscCall(MatSeqAIJGetArray(Y,&ya)); 3107 PetscCall(MatSeqAIJGetArrayRead(X,&xa)); 3108 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,xa,&one,ya,&one)); 3109 PetscCall(MatSeqAIJRestoreArrayRead(X,&xa)); 3110 PetscCall(MatSeqAIJRestoreArray(Y,&ya)); 3111 PetscCall(PetscLogFlops(2.0*bnz)); 3112 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 3113 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 3114 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 3115 PetscCall(MatAXPY_Basic(Y,a,X,str)); 3116 } else { 3117 Mat B; 3118 PetscInt *nnz; 3119 PetscCall(PetscMalloc1(Y->rmap->N,&nnz)); 3120 PetscCall(MatCreate(PetscObjectComm((PetscObject)Y),&B)); 3121 PetscCall(PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name)); 3122 PetscCall(MatSetLayouts(B,Y->rmap,Y->cmap)); 3123 PetscCall(MatSetType(B,((PetscObject)Y)->type_name)); 3124 PetscCall(MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz)); 3125 PetscCall(MatSeqAIJSetPreallocation(B,0,nnz)); 3126 PetscCall(MatAXPY_BasicWithPreallocation(B,Y,a,X,str)); 3127 PetscCall(MatHeaderMerge(Y,&B)); 3128 PetscCall(PetscFree(nnz)); 3129 } 3130 PetscFunctionReturn(0); 3131 } 3132 3133 PETSC_INTERN PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 3134 { 3135 #if defined(PETSC_USE_COMPLEX) 3136 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3137 PetscInt i,nz; 3138 PetscScalar *a; 3139 3140 PetscFunctionBegin; 3141 nz = aij->nz; 3142 PetscCall(MatSeqAIJGetArray(mat,&a)); 3143 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 3144 PetscCall(MatSeqAIJRestoreArray(mat,&a)); 3145 #else 3146 PetscFunctionBegin; 3147 #endif 3148 PetscFunctionReturn(0); 3149 } 3150 3151 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3152 { 3153 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3154 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3155 PetscReal atmp; 3156 PetscScalar *x; 3157 const MatScalar *aa,*av; 3158 3159 PetscFunctionBegin; 3160 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3161 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3162 aa = av; 3163 ai = a->i; 3164 aj = a->j; 3165 3166 PetscCall(VecSet(v,0.0)); 3167 PetscCall(VecGetArrayWrite(v,&x)); 3168 PetscCall(VecGetLocalSize(v,&n)); 3169 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3170 for (i=0; i<m; i++) { 3171 ncols = ai[1] - ai[0]; ai++; 3172 for (j=0; j<ncols; j++) { 3173 atmp = PetscAbsScalar(*aa); 3174 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3175 aa++; aj++; 3176 } 3177 } 3178 PetscCall(VecRestoreArrayWrite(v,&x)); 3179 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3180 PetscFunctionReturn(0); 3181 } 3182 3183 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3184 { 3185 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3186 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3187 PetscScalar *x; 3188 const MatScalar *aa,*av; 3189 3190 PetscFunctionBegin; 3191 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3192 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3193 aa = av; 3194 ai = a->i; 3195 aj = a->j; 3196 3197 PetscCall(VecSet(v,0.0)); 3198 PetscCall(VecGetArrayWrite(v,&x)); 3199 PetscCall(VecGetLocalSize(v,&n)); 3200 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3201 for (i=0; i<m; i++) { 3202 ncols = ai[1] - ai[0]; ai++; 3203 if (ncols == A->cmap->n) { /* row is dense */ 3204 x[i] = *aa; if (idx) idx[i] = 0; 3205 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3206 x[i] = 0.0; 3207 if (idx) { 3208 for (j=0; j<ncols; j++) { /* find first implicit 0.0 in the row */ 3209 if (aj[j] > j) { 3210 idx[i] = j; 3211 break; 3212 } 3213 } 3214 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3215 if (j==ncols && j < A->cmap->n) idx[i] = j; 3216 } 3217 } 3218 for (j=0; j<ncols; j++) { 3219 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3220 aa++; aj++; 3221 } 3222 } 3223 PetscCall(VecRestoreArrayWrite(v,&x)); 3224 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3225 PetscFunctionReturn(0); 3226 } 3227 3228 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3229 { 3230 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3231 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3232 PetscScalar *x; 3233 const MatScalar *aa,*av; 3234 3235 PetscFunctionBegin; 3236 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3237 aa = av; 3238 ai = a->i; 3239 aj = a->j; 3240 3241 PetscCall(VecSet(v,0.0)); 3242 PetscCall(VecGetArrayWrite(v,&x)); 3243 PetscCall(VecGetLocalSize(v,&n)); 3244 PetscCheck(n == m,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector, %" PetscInt_FMT " vs. %" PetscInt_FMT " rows", m, n); 3245 for (i=0; i<m; i++) { 3246 ncols = ai[1] - ai[0]; ai++; 3247 if (ncols == A->cmap->n) { /* row is dense */ 3248 x[i] = *aa; if (idx) idx[i] = 0; 3249 } else { /* row is sparse so already KNOW minimum is 0.0 or higher */ 3250 x[i] = 0.0; 3251 if (idx) { /* find first implicit 0.0 in the row */ 3252 for (j=0; j<ncols; j++) { 3253 if (aj[j] > j) { 3254 idx[i] = j; 3255 break; 3256 } 3257 } 3258 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3259 if (j==ncols && j < A->cmap->n) idx[i] = j; 3260 } 3261 } 3262 for (j=0; j<ncols; j++) { 3263 if (PetscAbsScalar(x[i]) > PetscAbsScalar(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3264 aa++; aj++; 3265 } 3266 } 3267 PetscCall(VecRestoreArrayWrite(v,&x)); 3268 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3269 PetscFunctionReturn(0); 3270 } 3271 3272 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3273 { 3274 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3275 PetscInt i,j,m = A->rmap->n,ncols,n; 3276 const PetscInt *ai,*aj; 3277 PetscScalar *x; 3278 const MatScalar *aa,*av; 3279 3280 PetscFunctionBegin; 3281 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3282 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3283 aa = av; 3284 ai = a->i; 3285 aj = a->j; 3286 3287 PetscCall(VecSet(v,0.0)); 3288 PetscCall(VecGetArrayWrite(v,&x)); 3289 PetscCall(VecGetLocalSize(v,&n)); 3290 PetscCheck(n == m,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3291 for (i=0; i<m; i++) { 3292 ncols = ai[1] - ai[0]; ai++; 3293 if (ncols == A->cmap->n) { /* row is dense */ 3294 x[i] = *aa; if (idx) idx[i] = 0; 3295 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3296 x[i] = 0.0; 3297 if (idx) { /* find first implicit 0.0 in the row */ 3298 for (j=0; j<ncols; j++) { 3299 if (aj[j] > j) { 3300 idx[i] = j; 3301 break; 3302 } 3303 } 3304 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3305 if (j==ncols && j < A->cmap->n) idx[i] = j; 3306 } 3307 } 3308 for (j=0; j<ncols; j++) { 3309 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3310 aa++; aj++; 3311 } 3312 } 3313 PetscCall(VecRestoreArrayWrite(v,&x)); 3314 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3315 PetscFunctionReturn(0); 3316 } 3317 3318 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3319 { 3320 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3321 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3322 MatScalar *diag,work[25],*v_work; 3323 const PetscReal shift = 0.0; 3324 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 3325 3326 PetscFunctionBegin; 3327 allowzeropivot = PetscNot(A->erroriffailure); 3328 if (a->ibdiagvalid) { 3329 if (values) *values = a->ibdiag; 3330 PetscFunctionReturn(0); 3331 } 3332 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 3333 if (!a->ibdiag) { 3334 PetscCall(PetscMalloc1(bs2*mbs,&a->ibdiag)); 3335 PetscCall(PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar))); 3336 } 3337 diag = a->ibdiag; 3338 if (values) *values = a->ibdiag; 3339 /* factor and invert each block */ 3340 switch (bs) { 3341 case 1: 3342 for (i=0; i<mbs; i++) { 3343 PetscCall(MatGetValues(A,1,&i,1,&i,diag+i)); 3344 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3345 if (allowzeropivot) { 3346 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3347 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3348 A->factorerror_zeropivot_row = i; 3349 PetscCall(PetscInfo(A,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON)); 3350 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON); 3351 } 3352 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3353 } 3354 break; 3355 case 2: 3356 for (i=0; i<mbs; i++) { 3357 ij[0] = 2*i; ij[1] = 2*i + 1; 3358 PetscCall(MatGetValues(A,2,ij,2,ij,diag)); 3359 PetscCall(PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected)); 3360 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3361 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 3362 diag += 4; 3363 } 3364 break; 3365 case 3: 3366 for (i=0; i<mbs; i++) { 3367 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3368 PetscCall(MatGetValues(A,3,ij,3,ij,diag)); 3369 PetscCall(PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected)); 3370 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3371 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 3372 diag += 9; 3373 } 3374 break; 3375 case 4: 3376 for (i=0; i<mbs; i++) { 3377 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3378 PetscCall(MatGetValues(A,4,ij,4,ij,diag)); 3379 PetscCall(PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected)); 3380 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3381 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 3382 diag += 16; 3383 } 3384 break; 3385 case 5: 3386 for (i=0; i<mbs; i++) { 3387 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3388 PetscCall(MatGetValues(A,5,ij,5,ij,diag)); 3389 PetscCall(PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 3390 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3391 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 3392 diag += 25; 3393 } 3394 break; 3395 case 6: 3396 for (i=0; i<mbs; i++) { 3397 ij[0] = 6*i; ij[1] = 6*i + 1; ij[2] = 6*i + 2; ij[3] = 6*i + 3; ij[4] = 6*i + 4; ij[5] = 6*i + 5; 3398 PetscCall(MatGetValues(A,6,ij,6,ij,diag)); 3399 PetscCall(PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected)); 3400 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3401 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 3402 diag += 36; 3403 } 3404 break; 3405 case 7: 3406 for (i=0; i<mbs; i++) { 3407 ij[0] = 7*i; ij[1] = 7*i + 1; ij[2] = 7*i + 2; ij[3] = 7*i + 3; ij[4] = 7*i + 4; ij[5] = 7*i + 5; ij[5] = 7*i + 6; 3408 PetscCall(MatGetValues(A,7,ij,7,ij,diag)); 3409 PetscCall(PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected)); 3410 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3411 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 3412 diag += 49; 3413 } 3414 break; 3415 default: 3416 PetscCall(PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ)); 3417 for (i=0; i<mbs; i++) { 3418 for (j=0; j<bs; j++) { 3419 IJ[j] = bs*i + j; 3420 } 3421 PetscCall(MatGetValues(A,bs,IJ,bs,IJ,diag)); 3422 PetscCall(PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 3423 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3424 PetscCall(PetscKernel_A_gets_transpose_A_N(diag,bs)); 3425 diag += bs2; 3426 } 3427 PetscCall(PetscFree3(v_work,v_pivots,IJ)); 3428 } 3429 a->ibdiagvalid = PETSC_TRUE; 3430 PetscFunctionReturn(0); 3431 } 3432 3433 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3434 { 3435 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3436 PetscScalar a,*aa; 3437 PetscInt m,n,i,j,col; 3438 3439 PetscFunctionBegin; 3440 if (!x->assembled) { 3441 PetscCall(MatGetSize(x,&m,&n)); 3442 for (i=0; i<m; i++) { 3443 for (j=0; j<aij->imax[i]; j++) { 3444 PetscCall(PetscRandomGetValue(rctx,&a)); 3445 col = (PetscInt)(n*PetscRealPart(a)); 3446 PetscCall(MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES)); 3447 } 3448 } 3449 } else { 3450 PetscCall(MatSeqAIJGetArrayWrite(x,&aa)); 3451 for (i=0; i<aij->nz; i++) PetscCall(PetscRandomGetValue(rctx,aa+i)); 3452 PetscCall(MatSeqAIJRestoreArrayWrite(x,&aa)); 3453 } 3454 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 3455 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 3456 PetscFunctionReturn(0); 3457 } 3458 3459 /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */ 3460 PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x,PetscInt low,PetscInt high,PetscRandom rctx) 3461 { 3462 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3463 PetscScalar a; 3464 PetscInt m,n,i,j,col,nskip; 3465 3466 PetscFunctionBegin; 3467 nskip = high - low; 3468 PetscCall(MatGetSize(x,&m,&n)); 3469 n -= nskip; /* shrink number of columns where nonzeros can be set */ 3470 for (i=0; i<m; i++) { 3471 for (j=0; j<aij->imax[i]; j++) { 3472 PetscCall(PetscRandomGetValue(rctx,&a)); 3473 col = (PetscInt)(n*PetscRealPart(a)); 3474 if (col >= low) col += nskip; /* shift col rightward to skip the hole */ 3475 PetscCall(MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES)); 3476 } 3477 } 3478 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 3479 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 3480 PetscFunctionReturn(0); 3481 } 3482 3483 /* -------------------------------------------------------------------*/ 3484 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3485 MatGetRow_SeqAIJ, 3486 MatRestoreRow_SeqAIJ, 3487 MatMult_SeqAIJ, 3488 /* 4*/ MatMultAdd_SeqAIJ, 3489 MatMultTranspose_SeqAIJ, 3490 MatMultTransposeAdd_SeqAIJ, 3491 NULL, 3492 NULL, 3493 NULL, 3494 /* 10*/ NULL, 3495 MatLUFactor_SeqAIJ, 3496 NULL, 3497 MatSOR_SeqAIJ, 3498 MatTranspose_SeqAIJ, 3499 /*1 5*/ MatGetInfo_SeqAIJ, 3500 MatEqual_SeqAIJ, 3501 MatGetDiagonal_SeqAIJ, 3502 MatDiagonalScale_SeqAIJ, 3503 MatNorm_SeqAIJ, 3504 /* 20*/ NULL, 3505 MatAssemblyEnd_SeqAIJ, 3506 MatSetOption_SeqAIJ, 3507 MatZeroEntries_SeqAIJ, 3508 /* 24*/ MatZeroRows_SeqAIJ, 3509 NULL, 3510 NULL, 3511 NULL, 3512 NULL, 3513 /* 29*/ MatSetUp_SeqAIJ, 3514 NULL, 3515 NULL, 3516 NULL, 3517 NULL, 3518 /* 34*/ MatDuplicate_SeqAIJ, 3519 NULL, 3520 NULL, 3521 MatILUFactor_SeqAIJ, 3522 NULL, 3523 /* 39*/ MatAXPY_SeqAIJ, 3524 MatCreateSubMatrices_SeqAIJ, 3525 MatIncreaseOverlap_SeqAIJ, 3526 MatGetValues_SeqAIJ, 3527 MatCopy_SeqAIJ, 3528 /* 44*/ MatGetRowMax_SeqAIJ, 3529 MatScale_SeqAIJ, 3530 MatShift_SeqAIJ, 3531 MatDiagonalSet_SeqAIJ, 3532 MatZeroRowsColumns_SeqAIJ, 3533 /* 49*/ MatSetRandom_SeqAIJ, 3534 MatGetRowIJ_SeqAIJ, 3535 MatRestoreRowIJ_SeqAIJ, 3536 MatGetColumnIJ_SeqAIJ, 3537 MatRestoreColumnIJ_SeqAIJ, 3538 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3539 NULL, 3540 NULL, 3541 MatPermute_SeqAIJ, 3542 NULL, 3543 /* 59*/ NULL, 3544 MatDestroy_SeqAIJ, 3545 MatView_SeqAIJ, 3546 NULL, 3547 NULL, 3548 /* 64*/ NULL, 3549 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3550 NULL, 3551 NULL, 3552 NULL, 3553 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3554 MatGetRowMinAbs_SeqAIJ, 3555 NULL, 3556 NULL, 3557 NULL, 3558 /* 74*/ NULL, 3559 MatFDColoringApply_AIJ, 3560 NULL, 3561 NULL, 3562 NULL, 3563 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3564 NULL, 3565 NULL, 3566 NULL, 3567 MatLoad_SeqAIJ, 3568 /* 84*/ MatIsSymmetric_SeqAIJ, 3569 MatIsHermitian_SeqAIJ, 3570 NULL, 3571 NULL, 3572 NULL, 3573 /* 89*/ NULL, 3574 NULL, 3575 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3576 NULL, 3577 NULL, 3578 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3579 NULL, 3580 NULL, 3581 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3582 NULL, 3583 /* 99*/ MatProductSetFromOptions_SeqAIJ, 3584 NULL, 3585 NULL, 3586 MatConjugate_SeqAIJ, 3587 NULL, 3588 /*104*/ MatSetValuesRow_SeqAIJ, 3589 MatRealPart_SeqAIJ, 3590 MatImaginaryPart_SeqAIJ, 3591 NULL, 3592 NULL, 3593 /*109*/ MatMatSolve_SeqAIJ, 3594 NULL, 3595 MatGetRowMin_SeqAIJ, 3596 NULL, 3597 MatMissingDiagonal_SeqAIJ, 3598 /*114*/ NULL, 3599 NULL, 3600 NULL, 3601 NULL, 3602 NULL, 3603 /*119*/ NULL, 3604 NULL, 3605 NULL, 3606 NULL, 3607 MatGetMultiProcBlock_SeqAIJ, 3608 /*124*/ MatFindNonzeroRows_SeqAIJ, 3609 MatGetColumnReductions_SeqAIJ, 3610 MatInvertBlockDiagonal_SeqAIJ, 3611 MatInvertVariableBlockDiagonal_SeqAIJ, 3612 NULL, 3613 /*129*/ NULL, 3614 NULL, 3615 NULL, 3616 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3617 MatTransposeColoringCreate_SeqAIJ, 3618 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3619 MatTransColoringApplyDenToSp_SeqAIJ, 3620 NULL, 3621 NULL, 3622 MatRARtNumeric_SeqAIJ_SeqAIJ, 3623 /*139*/NULL, 3624 NULL, 3625 NULL, 3626 MatFDColoringSetUp_SeqXAIJ, 3627 MatFindOffBlockDiagonalEntries_SeqAIJ, 3628 MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3629 /*145*/MatDestroySubMatrices_SeqAIJ, 3630 NULL, 3631 NULL 3632 }; 3633 3634 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3635 { 3636 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3637 PetscInt i,nz,n; 3638 3639 PetscFunctionBegin; 3640 nz = aij->maxnz; 3641 n = mat->rmap->n; 3642 for (i=0; i<nz; i++) { 3643 aij->j[i] = indices[i]; 3644 } 3645 aij->nz = nz; 3646 for (i=0; i<n; i++) { 3647 aij->ilen[i] = aij->imax[i]; 3648 } 3649 PetscFunctionReturn(0); 3650 } 3651 3652 /* 3653 * Given a sparse matrix with global column indices, compact it by using a local column space. 3654 * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3655 */ 3656 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3657 { 3658 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3659 PetscTable gid1_lid1; 3660 PetscTablePosition tpos; 3661 PetscInt gid,lid,i,ec,nz = aij->nz; 3662 PetscInt *garray,*jj = aij->j; 3663 3664 PetscFunctionBegin; 3665 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3666 PetscValidPointer(mapping,2); 3667 /* use a table */ 3668 PetscCall(PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1)); 3669 ec = 0; 3670 for (i=0; i<nz; i++) { 3671 PetscInt data,gid1 = jj[i] + 1; 3672 PetscCall(PetscTableFind(gid1_lid1,gid1,&data)); 3673 if (!data) { 3674 /* one based table */ 3675 PetscCall(PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES)); 3676 } 3677 } 3678 /* form array of columns we need */ 3679 PetscCall(PetscMalloc1(ec,&garray)); 3680 PetscCall(PetscTableGetHeadPosition(gid1_lid1,&tpos)); 3681 while (tpos) { 3682 PetscCall(PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid)); 3683 gid--; 3684 lid--; 3685 garray[lid] = gid; 3686 } 3687 PetscCall(PetscSortInt(ec,garray)); /* sort, and rebuild */ 3688 PetscCall(PetscTableRemoveAll(gid1_lid1)); 3689 for (i=0; i<ec; i++) { 3690 PetscCall(PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES)); 3691 } 3692 /* compact out the extra columns in B */ 3693 for (i=0; i<nz; i++) { 3694 PetscInt gid1 = jj[i] + 1; 3695 PetscCall(PetscTableFind(gid1_lid1,gid1,&lid)); 3696 lid--; 3697 jj[i] = lid; 3698 } 3699 PetscCall(PetscLayoutDestroy(&mat->cmap)); 3700 PetscCall(PetscTableDestroy(&gid1_lid1)); 3701 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat),ec,ec,1,&mat->cmap)); 3702 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping)); 3703 PetscCall(ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH)); 3704 PetscFunctionReturn(0); 3705 } 3706 3707 /*@ 3708 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3709 in the matrix. 3710 3711 Input Parameters: 3712 + mat - the SeqAIJ matrix 3713 - indices - the column indices 3714 3715 Level: advanced 3716 3717 Notes: 3718 This can be called if you have precomputed the nonzero structure of the 3719 matrix and want to provide it to the matrix object to improve the performance 3720 of the MatSetValues() operation. 3721 3722 You MUST have set the correct numbers of nonzeros per row in the call to 3723 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3724 3725 MUST be called before any calls to MatSetValues(); 3726 3727 The indices should start with zero, not one. 3728 3729 @*/ 3730 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3731 { 3732 PetscFunctionBegin; 3733 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3734 PetscValidIntPointer(indices,2); 3735 PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices)); 3736 PetscFunctionReturn(0); 3737 } 3738 3739 /* ----------------------------------------------------------------------------------------*/ 3740 3741 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3742 { 3743 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3744 size_t nz = aij->i[mat->rmap->n]; 3745 3746 PetscFunctionBegin; 3747 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3748 3749 /* allocate space for values if not already there */ 3750 if (!aij->saved_values) { 3751 PetscCall(PetscMalloc1(nz+1,&aij->saved_values)); 3752 PetscCall(PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar))); 3753 } 3754 3755 /* copy values over */ 3756 PetscCall(PetscArraycpy(aij->saved_values,aij->a,nz)); 3757 PetscFunctionReturn(0); 3758 } 3759 3760 /*@ 3761 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3762 example, reuse of the linear part of a Jacobian, while recomputing the 3763 nonlinear portion. 3764 3765 Collect on Mat 3766 3767 Input Parameters: 3768 . mat - the matrix (currently only AIJ matrices support this option) 3769 3770 Level: advanced 3771 3772 Common Usage, with SNESSolve(): 3773 $ Create Jacobian matrix 3774 $ Set linear terms into matrix 3775 $ Apply boundary conditions to matrix, at this time matrix must have 3776 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3777 $ boundary conditions again will not change the nonzero structure 3778 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3779 $ ierr = MatStoreValues(mat); 3780 $ Call SNESSetJacobian() with matrix 3781 $ In your Jacobian routine 3782 $ ierr = MatRetrieveValues(mat); 3783 $ Set nonlinear terms in matrix 3784 3785 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3786 $ // build linear portion of Jacobian 3787 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3788 $ ierr = MatStoreValues(mat); 3789 $ loop over nonlinear iterations 3790 $ ierr = MatRetrieveValues(mat); 3791 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3792 $ // call MatAssemblyBegin/End() on matrix 3793 $ Solve linear system with Jacobian 3794 $ endloop 3795 3796 Notes: 3797 Matrix must already be assemblied before calling this routine 3798 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3799 calling this routine. 3800 3801 When this is called multiple times it overwrites the previous set of stored values 3802 and does not allocated additional space. 3803 3804 .seealso: `MatRetrieveValues()` 3805 3806 @*/ 3807 PetscErrorCode MatStoreValues(Mat mat) 3808 { 3809 PetscFunctionBegin; 3810 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3811 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3812 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3813 PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat)); 3814 PetscFunctionReturn(0); 3815 } 3816 3817 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3818 { 3819 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3820 PetscInt nz = aij->i[mat->rmap->n]; 3821 3822 PetscFunctionBegin; 3823 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3824 PetscCheck(aij->saved_values,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3825 /* copy values over */ 3826 PetscCall(PetscArraycpy(aij->a,aij->saved_values,nz)); 3827 PetscFunctionReturn(0); 3828 } 3829 3830 /*@ 3831 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3832 example, reuse of the linear part of a Jacobian, while recomputing the 3833 nonlinear portion. 3834 3835 Collect on Mat 3836 3837 Input Parameters: 3838 . mat - the matrix (currently only AIJ matrices support this option) 3839 3840 Level: advanced 3841 3842 .seealso: `MatStoreValues()` 3843 3844 @*/ 3845 PetscErrorCode MatRetrieveValues(Mat mat) 3846 { 3847 PetscFunctionBegin; 3848 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3849 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3850 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3851 PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat)); 3852 PetscFunctionReturn(0); 3853 } 3854 3855 /* --------------------------------------------------------------------------------*/ 3856 /*@C 3857 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3858 (the default parallel PETSc format). For good matrix assembly performance 3859 the user should preallocate the matrix storage by setting the parameter nz 3860 (or the array nnz). By setting these parameters accurately, performance 3861 during matrix assembly can be increased by more than a factor of 50. 3862 3863 Collective 3864 3865 Input Parameters: 3866 + comm - MPI communicator, set to PETSC_COMM_SELF 3867 . m - number of rows 3868 . n - number of columns 3869 . nz - number of nonzeros per row (same for all rows) 3870 - nnz - array containing the number of nonzeros in the various rows 3871 (possibly different for each row) or NULL 3872 3873 Output Parameter: 3874 . A - the matrix 3875 3876 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3877 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3878 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3879 3880 Notes: 3881 If nnz is given then nz is ignored 3882 3883 The AIJ format (also called the Yale sparse matrix format or 3884 compressed row storage), is fully compatible with standard Fortran 77 3885 storage. That is, the stored row and column indices can begin at 3886 either one (as in Fortran) or zero. See the users' manual for details. 3887 3888 Specify the preallocated storage with either nz or nnz (not both). 3889 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3890 allocation. For large problems you MUST preallocate memory or you 3891 will get TERRIBLE performance, see the users' manual chapter on matrices. 3892 3893 By default, this format uses inodes (identical nodes) when possible, to 3894 improve numerical efficiency of matrix-vector products and solves. We 3895 search for consecutive rows with the same nonzero structure, thereby 3896 reusing matrix information to achieve increased efficiency. 3897 3898 Options Database Keys: 3899 + -mat_no_inode - Do not use inodes 3900 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3901 3902 Level: intermediate 3903 3904 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()` 3905 3906 @*/ 3907 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3908 { 3909 PetscFunctionBegin; 3910 PetscCall(MatCreate(comm,A)); 3911 PetscCall(MatSetSizes(*A,m,n,m,n)); 3912 PetscCall(MatSetType(*A,MATSEQAIJ)); 3913 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz)); 3914 PetscFunctionReturn(0); 3915 } 3916 3917 /*@C 3918 MatSeqAIJSetPreallocation - For good matrix assembly performance 3919 the user should preallocate the matrix storage by setting the parameter nz 3920 (or the array nnz). By setting these parameters accurately, performance 3921 during matrix assembly can be increased by more than a factor of 50. 3922 3923 Collective 3924 3925 Input Parameters: 3926 + B - The matrix 3927 . nz - number of nonzeros per row (same for all rows) 3928 - nnz - array containing the number of nonzeros in the various rows 3929 (possibly different for each row) or NULL 3930 3931 Notes: 3932 If nnz is given then nz is ignored 3933 3934 The AIJ format (also called the Yale sparse matrix format or 3935 compressed row storage), is fully compatible with standard Fortran 77 3936 storage. That is, the stored row and column indices can begin at 3937 either one (as in Fortran) or zero. See the users' manual for details. 3938 3939 Specify the preallocated storage with either nz or nnz (not both). 3940 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3941 allocation. For large problems you MUST preallocate memory or you 3942 will get TERRIBLE performance, see the users' manual chapter on matrices. 3943 3944 You can call MatGetInfo() to get information on how effective the preallocation was; 3945 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3946 You can also run with the option -info and look for messages with the string 3947 malloc in them to see if additional memory allocation was needed. 3948 3949 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3950 entries or columns indices 3951 3952 By default, this format uses inodes (identical nodes) when possible, to 3953 improve numerical efficiency of matrix-vector products and solves. We 3954 search for consecutive rows with the same nonzero structure, thereby 3955 reusing matrix information to achieve increased efficiency. 3956 3957 Options Database Keys: 3958 + -mat_no_inode - Do not use inodes 3959 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3960 3961 Level: intermediate 3962 3963 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()`, `MatGetInfo()`, 3964 `MatSeqAIJSetTotalPreallocation()` 3965 3966 @*/ 3967 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3968 { 3969 PetscFunctionBegin; 3970 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3971 PetscValidType(B,1); 3972 PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz)); 3973 PetscFunctionReturn(0); 3974 } 3975 3976 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3977 { 3978 Mat_SeqAIJ *b; 3979 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3980 PetscInt i; 3981 3982 PetscFunctionBegin; 3983 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3984 if (nz == MAT_SKIP_ALLOCATION) { 3985 skipallocation = PETSC_TRUE; 3986 nz = 0; 3987 } 3988 PetscCall(PetscLayoutSetUp(B->rmap)); 3989 PetscCall(PetscLayoutSetUp(B->cmap)); 3990 3991 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3992 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %" PetscInt_FMT,nz); 3993 if (PetscUnlikelyDebug(nnz)) { 3994 for (i=0; i<B->rmap->n; i++) { 3995 PetscCheck(nnz[i] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %" PetscInt_FMT " value %" PetscInt_FMT,i,nnz[i]); 3996 PetscCheck(nnz[i] <= B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " rowlength %" PetscInt_FMT,i,nnz[i],B->cmap->n); 3997 } 3998 } 3999 4000 B->preallocated = PETSC_TRUE; 4001 4002 b = (Mat_SeqAIJ*)B->data; 4003 4004 if (!skipallocation) { 4005 if (!b->imax) { 4006 PetscCall(PetscMalloc1(B->rmap->n,&b->imax)); 4007 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4008 } 4009 if (!b->ilen) { 4010 /* b->ilen will count nonzeros in each row so far. */ 4011 PetscCall(PetscCalloc1(B->rmap->n,&b->ilen)); 4012 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4013 } else { 4014 PetscCall(PetscMemzero(b->ilen,B->rmap->n*sizeof(PetscInt))); 4015 } 4016 if (!b->ipre) { 4017 PetscCall(PetscMalloc1(B->rmap->n,&b->ipre)); 4018 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4019 } 4020 if (!nnz) { 4021 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 4022 else if (nz < 0) nz = 1; 4023 nz = PetscMin(nz,B->cmap->n); 4024 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 4025 nz = nz*B->rmap->n; 4026 } else { 4027 PetscInt64 nz64 = 0; 4028 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 4029 PetscCall(PetscIntCast(nz64,&nz)); 4030 } 4031 4032 /* allocate the matrix space */ 4033 /* FIXME: should B's old memory be unlogged? */ 4034 PetscCall(MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i)); 4035 if (B->structure_only) { 4036 PetscCall(PetscMalloc1(nz,&b->j)); 4037 PetscCall(PetscMalloc1(B->rmap->n+1,&b->i)); 4038 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt))); 4039 } else { 4040 PetscCall(PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i)); 4041 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)))); 4042 } 4043 b->i[0] = 0; 4044 for (i=1; i<B->rmap->n+1; i++) { 4045 b->i[i] = b->i[i-1] + b->imax[i-1]; 4046 } 4047 if (B->structure_only) { 4048 b->singlemalloc = PETSC_FALSE; 4049 b->free_a = PETSC_FALSE; 4050 } else { 4051 b->singlemalloc = PETSC_TRUE; 4052 b->free_a = PETSC_TRUE; 4053 } 4054 b->free_ij = PETSC_TRUE; 4055 } else { 4056 b->free_a = PETSC_FALSE; 4057 b->free_ij = PETSC_FALSE; 4058 } 4059 4060 if (b->ipre && nnz != b->ipre && b->imax) { 4061 /* reserve user-requested sparsity */ 4062 PetscCall(PetscArraycpy(b->ipre,b->imax,B->rmap->n)); 4063 } 4064 4065 b->nz = 0; 4066 b->maxnz = nz; 4067 B->info.nz_unneeded = (double)b->maxnz; 4068 if (realalloc) PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE)); 4069 B->was_assembled = PETSC_FALSE; 4070 B->assembled = PETSC_FALSE; 4071 /* We simply deem preallocation has changed nonzero state. Updating the state 4072 will give clients (like AIJKokkos) a chance to know something has happened. 4073 */ 4074 B->nonzerostate++; 4075 PetscFunctionReturn(0); 4076 } 4077 4078 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 4079 { 4080 Mat_SeqAIJ *a; 4081 PetscInt i; 4082 4083 PetscFunctionBegin; 4084 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4085 4086 /* Check local size. If zero, then return */ 4087 if (!A->rmap->n) PetscFunctionReturn(0); 4088 4089 a = (Mat_SeqAIJ*)A->data; 4090 /* if no saved info, we error out */ 4091 PetscCheck(a->ipre,PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"No saved preallocation info "); 4092 4093 PetscCheck(a->i && a->j && a->a && a->imax && a->ilen,PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Memory info is incomplete, and can not reset preallocation "); 4094 4095 PetscCall(PetscArraycpy(a->imax,a->ipre,A->rmap->n)); 4096 PetscCall(PetscArrayzero(a->ilen,A->rmap->n)); 4097 a->i[0] = 0; 4098 for (i=1; i<A->rmap->n+1; i++) { 4099 a->i[i] = a->i[i-1] + a->imax[i-1]; 4100 } 4101 A->preallocated = PETSC_TRUE; 4102 a->nz = 0; 4103 a->maxnz = a->i[A->rmap->n]; 4104 A->info.nz_unneeded = (double)a->maxnz; 4105 A->was_assembled = PETSC_FALSE; 4106 A->assembled = PETSC_FALSE; 4107 PetscFunctionReturn(0); 4108 } 4109 4110 /*@ 4111 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 4112 4113 Input Parameters: 4114 + B - the matrix 4115 . i - the indices into j for the start of each row (starts with zero) 4116 . j - the column indices for each row (starts with zero) these must be sorted for each row 4117 - v - optional values in the matrix 4118 4119 Level: developer 4120 4121 Notes: 4122 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 4123 4124 This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero 4125 structure will be the union of all the previous nonzero structures. 4126 4127 Developer Notes: 4128 An optimization could be added to the implementation where it checks if the i, and j are identical to the current i and j and 4129 then just copies the v values directly with PetscMemcpy(). 4130 4131 This routine could also take a PetscCopyMode argument to allow sharing the values instead of always copying them. 4132 4133 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatSeqAIJSetPreallocation()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MatResetPreallocation()` 4134 @*/ 4135 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 4136 { 4137 PetscFunctionBegin; 4138 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 4139 PetscValidType(B,1); 4140 PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v)); 4141 PetscFunctionReturn(0); 4142 } 4143 4144 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 4145 { 4146 PetscInt i; 4147 PetscInt m,n; 4148 PetscInt nz; 4149 PetscInt *nnz; 4150 4151 PetscFunctionBegin; 4152 PetscCheck(Ii[0] == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %" PetscInt_FMT, Ii[0]); 4153 4154 PetscCall(PetscLayoutSetUp(B->rmap)); 4155 PetscCall(PetscLayoutSetUp(B->cmap)); 4156 4157 PetscCall(MatGetSize(B, &m, &n)); 4158 PetscCall(PetscMalloc1(m+1, &nnz)); 4159 for (i = 0; i < m; i++) { 4160 nz = Ii[i+1]- Ii[i]; 4161 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 4162 nnz[i] = nz; 4163 } 4164 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 4165 PetscCall(PetscFree(nnz)); 4166 4167 for (i = 0; i < m; i++) { 4168 PetscCall(MatSetValues_SeqAIJ(B, 1, &i, Ii[i+1] - Ii[i], J+Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES)); 4169 } 4170 4171 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 4172 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 4173 4174 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 4175 PetscFunctionReturn(0); 4176 } 4177 4178 /*@ 4179 MatSeqAIJKron - Computes C, the Kronecker product of A and B. 4180 4181 Input Parameters: 4182 + A - left-hand side matrix 4183 . B - right-hand side matrix 4184 - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 4185 4186 Output Parameter: 4187 . C - Kronecker product of A and B 4188 4189 Level: intermediate 4190 4191 Notes: 4192 MAT_REUSE_MATRIX can only be used when the nonzero structure of the product matrix has not changed from that last call to MatSeqAIJKron(). 4193 4194 .seealso: `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATKAIJ`, `MatReuse` 4195 @*/ 4196 PetscErrorCode MatSeqAIJKron(Mat A,Mat B,MatReuse reuse,Mat *C) 4197 { 4198 PetscFunctionBegin; 4199 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4200 PetscValidType(A,1); 4201 PetscValidHeaderSpecific(B,MAT_CLASSID,2); 4202 PetscValidType(B,2); 4203 PetscValidPointer(C,4); 4204 if (reuse == MAT_REUSE_MATRIX) { 4205 PetscValidHeaderSpecific(*C,MAT_CLASSID,4); 4206 PetscValidType(*C,4); 4207 } 4208 PetscTryMethod(A,"MatSeqAIJKron_C",(Mat,Mat,MatReuse,Mat*),(A,B,reuse,C)); 4209 PetscFunctionReturn(0); 4210 } 4211 4212 PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A,Mat B,MatReuse reuse,Mat *C) 4213 { 4214 Mat newmat; 4215 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4216 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 4217 PetscScalar *v; 4218 const PetscScalar *aa,*ba; 4219 PetscInt *i,*j,m,n,p,q,nnz = 0,am = A->rmap->n,bm = B->rmap->n,an = A->cmap->n, bn = B->cmap->n; 4220 PetscBool flg; 4221 4222 PetscFunctionBegin; 4223 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4224 PetscCheck(A->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4225 PetscCheck(!B->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4226 PetscCheck(B->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4227 PetscCall(PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&flg)); 4228 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatType %s",((PetscObject)B)->type_name); 4229 PetscCheck(reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatReuse %d",(int)reuse); 4230 if (reuse == MAT_INITIAL_MATRIX) { 4231 PetscCall(PetscMalloc2(am*bm+1,&i,a->i[am]*b->i[bm],&j)); 4232 PetscCall(MatCreate(PETSC_COMM_SELF,&newmat)); 4233 PetscCall(MatSetSizes(newmat,am*bm,an*bn,am*bm,an*bn)); 4234 PetscCall(MatSetType(newmat,MATAIJ)); 4235 i[0] = 0; 4236 for (m = 0; m < am; ++m) { 4237 for (p = 0; p < bm; ++p) { 4238 i[m*bm + p + 1] = i[m*bm + p] + (a->i[m+1] - a->i[m]) * (b->i[p+1] - b->i[p]); 4239 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4240 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4241 j[nnz++] = a->j[n]*bn + b->j[q]; 4242 } 4243 } 4244 } 4245 } 4246 PetscCall(MatSeqAIJSetPreallocationCSR(newmat,i,j,NULL)); 4247 *C = newmat; 4248 PetscCall(PetscFree2(i,j)); 4249 nnz = 0; 4250 } 4251 PetscCall(MatSeqAIJGetArray(*C,&v)); 4252 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4253 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 4254 for (m = 0; m < am; ++m) { 4255 for (p = 0; p < bm; ++p) { 4256 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4257 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4258 v[nnz++] = aa[n] * ba[q]; 4259 } 4260 } 4261 } 4262 } 4263 PetscCall(MatSeqAIJRestoreArray(*C,&v)); 4264 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 4265 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 4266 PetscFunctionReturn(0); 4267 } 4268 4269 #include <../src/mat/impls/dense/seq/dense.h> 4270 #include <petsc/private/kernels/petscaxpy.h> 4271 4272 /* 4273 Computes (B'*A')' since computing B*A directly is untenable 4274 4275 n p p 4276 [ ] [ ] [ ] 4277 m [ A ] * n [ B ] = m [ C ] 4278 [ ] [ ] [ ] 4279 4280 */ 4281 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 4282 { 4283 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 4284 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 4285 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 4286 PetscInt i,j,n,m,q,p; 4287 const PetscInt *ii,*idx; 4288 const PetscScalar *b,*a,*a_q; 4289 PetscScalar *c,*c_q; 4290 PetscInt clda = sub_c->lda; 4291 PetscInt alda = sub_a->lda; 4292 4293 PetscFunctionBegin; 4294 m = A->rmap->n; 4295 n = A->cmap->n; 4296 p = B->cmap->n; 4297 a = sub_a->v; 4298 b = sub_b->a; 4299 c = sub_c->v; 4300 if (clda == m) { 4301 PetscCall(PetscArrayzero(c,m*p)); 4302 } else { 4303 for (j=0;j<p;j++) 4304 for (i=0;i<m;i++) 4305 c[j*clda + i] = 0.0; 4306 } 4307 ii = sub_b->i; 4308 idx = sub_b->j; 4309 for (i=0; i<n; i++) { 4310 q = ii[i+1] - ii[i]; 4311 while (q-->0) { 4312 c_q = c + clda*(*idx); 4313 a_q = a + alda*i; 4314 PetscKernelAXPY(c_q,*b,a_q,m); 4315 idx++; 4316 b++; 4317 } 4318 } 4319 PetscFunctionReturn(0); 4320 } 4321 4322 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C) 4323 { 4324 PetscInt m=A->rmap->n,n=B->cmap->n; 4325 PetscBool cisdense; 4326 4327 PetscFunctionBegin; 4328 PetscCheck(A->cmap->n == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %" PetscInt_FMT " != B->rmap->n %" PetscInt_FMT,A->cmap->n,B->rmap->n); 4329 PetscCall(MatSetSizes(C,m,n,m,n)); 4330 PetscCall(MatSetBlockSizesFromMats(C,A,B)); 4331 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 4332 if (!cisdense) { 4333 PetscCall(MatSetType(C,MATDENSE)); 4334 } 4335 PetscCall(MatSetUp(C)); 4336 4337 C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4338 PetscFunctionReturn(0); 4339 } 4340 4341 /* ----------------------------------------------------------------*/ 4342 /*MC 4343 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4344 based on compressed sparse row format. 4345 4346 Options Database Keys: 4347 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4348 4349 Level: beginner 4350 4351 Notes: 4352 MatSetValues() may be called for this matrix type with a NULL argument for the numerical values, 4353 in this case the values associated with the rows and columns one passes in are set to zero 4354 in the matrix 4355 4356 MatSetOptions(,MAT_STRUCTURE_ONLY,PETSC_TRUE) may be called for this matrix type. In this no 4357 space is allocated for the nonzero entries and any entries passed with MatSetValues() are ignored 4358 4359 Developer Notes: 4360 It would be nice if all matrix formats supported passing NULL in for the numerical values 4361 4362 .seealso: `MatCreateSeqAIJ()`, `MatSetFromOptions()`, `MatSetType()`, `MatCreate()`, `MatType`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4363 M*/ 4364 4365 /*MC 4366 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4367 4368 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 4369 and MATMPIAIJ otherwise. As a result, for single process communicators, 4370 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4371 for communicators controlling multiple processes. It is recommended that you call both of 4372 the above preallocation routines for simplicity. 4373 4374 Options Database Keys: 4375 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 4376 4377 Developer Notes: 4378 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 4379 enough exist. 4380 4381 Level: beginner 4382 4383 .seealso: `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATMPIAIJ`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4384 M*/ 4385 4386 /*MC 4387 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4388 4389 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 4390 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 4391 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4392 for communicators controlling multiple processes. It is recommended that you call both of 4393 the above preallocation routines for simplicity. 4394 4395 Options Database Keys: 4396 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 4397 4398 Level: beginner 4399 4400 .seealso: `MatCreateMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL` 4401 M*/ 4402 4403 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 4404 #if defined(PETSC_HAVE_ELEMENTAL) 4405 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 4406 #endif 4407 #if defined(PETSC_HAVE_SCALAPACK) 4408 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*); 4409 #endif 4410 #if defined(PETSC_HAVE_HYPRE) 4411 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 4412 #endif 4413 4414 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 4415 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 4416 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat); 4417 4418 /*@C 4419 MatSeqAIJGetArray - gives read/write access to the array where the data for a MATSEQAIJ matrix is stored 4420 4421 Not Collective 4422 4423 Input Parameter: 4424 . mat - a MATSEQAIJ matrix 4425 4426 Output Parameter: 4427 . array - pointer to the data 4428 4429 Level: intermediate 4430 4431 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4432 @*/ 4433 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4434 { 4435 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4436 4437 PetscFunctionBegin; 4438 if (aij->ops->getarray) { 4439 PetscCall((*aij->ops->getarray)(A,array)); 4440 } else { 4441 *array = aij->a; 4442 } 4443 PetscFunctionReturn(0); 4444 } 4445 4446 /*@C 4447 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4448 4449 Not Collective 4450 4451 Input Parameters: 4452 + mat - a MATSEQAIJ matrix 4453 - array - pointer to the data 4454 4455 Level: intermediate 4456 4457 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayF90()` 4458 @*/ 4459 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4460 { 4461 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4462 4463 PetscFunctionBegin; 4464 if (aij->ops->restorearray) { 4465 PetscCall((*aij->ops->restorearray)(A,array)); 4466 } else { 4467 *array = NULL; 4468 } 4469 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4470 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4471 PetscFunctionReturn(0); 4472 } 4473 4474 /*@C 4475 MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a MATSEQAIJ matrix is stored 4476 4477 Not Collective 4478 4479 Input Parameter: 4480 . mat - a MATSEQAIJ matrix 4481 4482 Output Parameter: 4483 . array - pointer to the data 4484 4485 Level: intermediate 4486 4487 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4488 @*/ 4489 PetscErrorCode MatSeqAIJGetArrayRead(Mat A,const PetscScalar **array) 4490 { 4491 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4492 4493 PetscFunctionBegin; 4494 if (aij->ops->getarrayread) { 4495 PetscCall((*aij->ops->getarrayread)(A,array)); 4496 } else { 4497 *array = aij->a; 4498 } 4499 PetscFunctionReturn(0); 4500 } 4501 4502 /*@C 4503 MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4504 4505 Not Collective 4506 4507 Input Parameter: 4508 . mat - a MATSEQAIJ matrix 4509 4510 Output Parameter: 4511 . array - pointer to the data 4512 4513 Level: intermediate 4514 4515 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4516 @*/ 4517 PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A,const PetscScalar **array) 4518 { 4519 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4520 4521 PetscFunctionBegin; 4522 if (aij->ops->restorearrayread) { 4523 PetscCall((*aij->ops->restorearrayread)(A,array)); 4524 } else { 4525 *array = NULL; 4526 } 4527 PetscFunctionReturn(0); 4528 } 4529 4530 /*@C 4531 MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a MATSEQAIJ matrix is stored 4532 4533 Not Collective 4534 4535 Input Parameter: 4536 . mat - a MATSEQAIJ matrix 4537 4538 Output Parameter: 4539 . array - pointer to the data 4540 4541 Level: intermediate 4542 4543 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4544 @*/ 4545 PetscErrorCode MatSeqAIJGetArrayWrite(Mat A,PetscScalar **array) 4546 { 4547 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4548 4549 PetscFunctionBegin; 4550 if (aij->ops->getarraywrite) { 4551 PetscCall((*aij->ops->getarraywrite)(A,array)); 4552 } else { 4553 *array = aij->a; 4554 } 4555 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4556 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4557 PetscFunctionReturn(0); 4558 } 4559 4560 /*@C 4561 MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4562 4563 Not Collective 4564 4565 Input Parameter: 4566 . mat - a MATSEQAIJ matrix 4567 4568 Output Parameter: 4569 . array - pointer to the data 4570 4571 Level: intermediate 4572 4573 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4574 @*/ 4575 PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A,PetscScalar **array) 4576 { 4577 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4578 4579 PetscFunctionBegin; 4580 if (aij->ops->restorearraywrite) { 4581 PetscCall((*aij->ops->restorearraywrite)(A,array)); 4582 } else { 4583 *array = NULL; 4584 } 4585 PetscFunctionReturn(0); 4586 } 4587 4588 /*@C 4589 MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the SEQAIJ matrix 4590 4591 Not Collective 4592 4593 Input Parameter: 4594 . mat - a matrix of type MATSEQAIJ or its subclasses 4595 4596 Output Parameters: 4597 + i - row map array of the matrix 4598 . j - column index array of the matrix 4599 . a - data array of the matrix 4600 - memtype - memory type of the arrays 4601 4602 Notes: 4603 Any of the output parameters can be NULL, in which case the corresponding value is not returned. 4604 If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host. 4605 4606 One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix. 4607 If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix. 4608 4609 Level: Developer 4610 4611 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4612 @*/ 4613 PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat,const PetscInt **i,const PetscInt **j,PetscScalar **a,PetscMemType *mtype) 4614 { 4615 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 4616 4617 PetscFunctionBegin; 4618 PetscCheck(mat->preallocated,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"matrix is not preallocated"); 4619 if (aij->ops->getcsrandmemtype) { 4620 PetscCall((*aij->ops->getcsrandmemtype)(mat,i,j,a,mtype)); 4621 } else { 4622 if (i) *i = aij->i; 4623 if (j) *j = aij->j; 4624 if (a) *a = aij->a; 4625 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 4626 } 4627 PetscFunctionReturn(0); 4628 } 4629 4630 /*@C 4631 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4632 4633 Not Collective 4634 4635 Input Parameter: 4636 . mat - a MATSEQAIJ matrix 4637 4638 Output Parameter: 4639 . nz - the maximum number of nonzeros in any row 4640 4641 Level: intermediate 4642 4643 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4644 @*/ 4645 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4646 { 4647 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4648 4649 PetscFunctionBegin; 4650 *nz = aij->rmax; 4651 PetscFunctionReturn(0); 4652 } 4653 4654 PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, const PetscInt coo_i[], const PetscInt coo_j[]) 4655 { 4656 MPI_Comm comm; 4657 PetscInt *i,*j; 4658 PetscInt M,N,row; 4659 PetscCount k,p,q,nneg,nnz,start,end; /* Index the coo array, so use PetscCount as their type */ 4660 PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */ 4661 PetscInt *Aj; 4662 PetscScalar *Aa; 4663 Mat_SeqAIJ *seqaij = (Mat_SeqAIJ*)(mat->data); 4664 MatType rtype; 4665 PetscCount *perm,*jmap; 4666 4667 PetscFunctionBegin; 4668 PetscCall(MatResetPreallocationCOO_SeqAIJ(mat)); 4669 PetscCall(PetscObjectGetComm((PetscObject)mat,&comm)); 4670 PetscCall(MatGetSize(mat,&M,&N)); 4671 PetscCall(PetscMalloc2(coo_n,&i,coo_n,&j)); 4672 PetscCall(PetscArraycpy(i,coo_i,coo_n)); /* Make a copy since we'll modify it */ 4673 PetscCall(PetscArraycpy(j,coo_j,coo_n)); 4674 PetscCall(PetscMalloc1(coo_n,&perm)); 4675 for (k=0; k<coo_n; k++) { /* Ignore entries with negative row or col indices */ 4676 if (j[k] < 0) i[k] = -1; 4677 perm[k] = k; 4678 } 4679 4680 /* Sort by row */ 4681 PetscCall(PetscSortIntWithIntCountArrayPair(coo_n,i,j,perm)); 4682 for (k=0; k<coo_n; k++) {if (i[k] >= 0) break;} /* Advance k to the first row with a non-negative index */ 4683 nneg = k; 4684 PetscCall(PetscMalloc1(coo_n-nneg+1,&jmap)); /* +1 to make a CSR-like data structure. jmap[i] originally is the number of repeats for i-th nonzero */ 4685 nnz = 0; /* Total number of unique nonzeros to be counted */ 4686 jmap++; /* Inc jmap by 1 for convinience */ 4687 4688 PetscCall(PetscCalloc1(M+1,&Ai)); /* CSR of A */ 4689 PetscCall(PetscMalloc1(coo_n-nneg,&Aj)); /* We have at most coo_n-nneg unique nonzeros */ 4690 4691 /* In each row, sort by column, then unique column indices to get row length */ 4692 Ai++; /* Inc by 1 for convinience */ 4693 q = 0; /* q-th unique nonzero, with q starting from 0 */ 4694 while (k<coo_n) { 4695 row = i[k]; 4696 start = k; /* [start,end) indices for this row */ 4697 while (k<coo_n && i[k] == row) k++; 4698 end = k; 4699 PetscCall(PetscSortIntWithCountArray(end-start,j+start,perm+start)); 4700 /* Find number of unique col entries in this row */ 4701 Aj[q] = j[start]; /* Log the first nonzero in this row */ 4702 jmap[q] = 1; /* Number of repeats of this nozero entry */ 4703 Ai[row] = 1; 4704 nnz++; 4705 4706 for (p=start+1; p<end; p++) { /* Scan remaining nonzero in this row */ 4707 if (j[p] != j[p-1]) { /* Meet a new nonzero */ 4708 q++; 4709 jmap[q] = 1; 4710 Aj[q] = j[p]; 4711 Ai[row]++; 4712 nnz++; 4713 } else { 4714 jmap[q]++; 4715 } 4716 } 4717 q++; /* Move to next row and thus next unique nonzero */ 4718 } 4719 PetscCall(PetscFree2(i,j)); 4720 4721 Ai--; /* Back to the beginning of Ai[] */ 4722 for (k=0; k<M; k++) Ai[k+1] += Ai[k]; 4723 jmap--; /* Back to the beginning of jmap[] */ 4724 jmap[0] = 0; 4725 for (k=0; k<nnz; k++) jmap[k+1] += jmap[k]; 4726 if (nnz < coo_n-nneg) { /* Realloc with actual number of unique nonzeros */ 4727 PetscCount *jmap_new; 4728 PetscInt *Aj_new; 4729 4730 PetscCall(PetscMalloc1(nnz+1,&jmap_new)); 4731 PetscCall(PetscArraycpy(jmap_new,jmap,nnz+1)); 4732 PetscCall(PetscFree(jmap)); 4733 jmap = jmap_new; 4734 4735 PetscCall(PetscMalloc1(nnz,&Aj_new)); 4736 PetscCall(PetscArraycpy(Aj_new,Aj,nnz)); 4737 PetscCall(PetscFree(Aj)); 4738 Aj = Aj_new; 4739 } 4740 4741 if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */ 4742 PetscCount *perm_new; 4743 4744 PetscCall(PetscMalloc1(coo_n-nneg,&perm_new)); 4745 PetscCall(PetscArraycpy(perm_new,perm+nneg,coo_n-nneg)); 4746 PetscCall(PetscFree(perm)); 4747 perm = perm_new; 4748 } 4749 4750 PetscCall(MatGetRootType_Private(mat,&rtype)); 4751 PetscCall(PetscCalloc1(nnz,&Aa)); /* Zero the matrix */ 4752 PetscCall(MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF,M,N,Ai,Aj,Aa,rtype,mat)); 4753 4754 seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */ 4755 seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */ 4756 /* Record COO fields */ 4757 seqaij->coo_n = coo_n; 4758 seqaij->Atot = coo_n-nneg; /* Annz is seqaij->nz, so no need to record that again */ 4759 seqaij->jmap = jmap; /* of length nnz+1 */ 4760 seqaij->perm = perm; 4761 PetscFunctionReturn(0); 4762 } 4763 4764 static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A,const PetscScalar v[],InsertMode imode) 4765 { 4766 Mat_SeqAIJ *aseq = (Mat_SeqAIJ*)A->data; 4767 PetscCount i,j,Annz = aseq->nz; 4768 PetscCount *perm = aseq->perm,*jmap = aseq->jmap; 4769 PetscScalar *Aa; 4770 4771 PetscFunctionBegin; 4772 PetscCall(MatSeqAIJGetArray(A,&Aa)); 4773 for (i=0; i<Annz; i++) { 4774 PetscScalar sum = 0.0; 4775 for (j=jmap[i]; j<jmap[i+1]; j++) sum += v[perm[j]]; 4776 Aa[i] = (imode == INSERT_VALUES? 0.0 : Aa[i]) + sum; 4777 } 4778 PetscCall(MatSeqAIJRestoreArray(A,&Aa)); 4779 PetscFunctionReturn(0); 4780 } 4781 4782 #if defined(PETSC_HAVE_CUDA) 4783 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat,MatType,MatReuse,Mat*); 4784 #endif 4785 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4786 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat,MatType,MatReuse,Mat*); 4787 #endif 4788 4789 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4790 { 4791 Mat_SeqAIJ *b; 4792 PetscMPIInt size; 4793 4794 PetscFunctionBegin; 4795 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 4796 PetscCheck(size <= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4797 4798 PetscCall(PetscNewLog(B,&b)); 4799 4800 B->data = (void*)b; 4801 4802 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 4803 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4804 4805 b->row = NULL; 4806 b->col = NULL; 4807 b->icol = NULL; 4808 b->reallocs = 0; 4809 b->ignorezeroentries = PETSC_FALSE; 4810 b->roworiented = PETSC_TRUE; 4811 b->nonew = 0; 4812 b->diag = NULL; 4813 b->solve_work = NULL; 4814 B->spptr = NULL; 4815 b->saved_values = NULL; 4816 b->idiag = NULL; 4817 b->mdiag = NULL; 4818 b->ssor_work = NULL; 4819 b->omega = 1.0; 4820 b->fshift = 0.0; 4821 b->idiagvalid = PETSC_FALSE; 4822 b->ibdiagvalid = PETSC_FALSE; 4823 b->keepnonzeropattern = PETSC_FALSE; 4824 4825 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4826 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4827 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ)); 4828 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ)); 4829 #endif 4830 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ)); 4831 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ)); 4832 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ)); 4833 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ)); 4834 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ)); 4835 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM)); 4836 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL)); 4837 #if defined(PETSC_HAVE_MKL_SPARSE) 4838 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL)); 4839 #endif 4840 #if defined(PETSC_HAVE_CUDA) 4841 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 4842 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4843 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",MatProductSetFromOptions_SeqAIJ)); 4844 #endif 4845 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4846 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijkokkos_C",MatConvert_SeqAIJ_SeqAIJKokkos)); 4847 #endif 4848 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL)); 4849 #if defined(PETSC_HAVE_ELEMENTAL) 4850 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental)); 4851 #endif 4852 #if defined(PETSC_HAVE_SCALAPACK) 4853 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_scalapack_C",MatConvert_AIJ_ScaLAPACK)); 4854 #endif 4855 #if defined(PETSC_HAVE_HYPRE) 4856 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE)); 4857 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",MatProductSetFromOptions_Transpose_AIJ_AIJ)); 4858 #endif 4859 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense)); 4860 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL)); 4861 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS)); 4862 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ)); 4863 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ)); 4864 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ)); 4865 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ)); 4866 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ)); 4867 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ)); 4868 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_is_seqaij_C",MatProductSetFromOptions_IS_XAIJ)); 4869 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqaij_C",MatProductSetFromOptions_SeqDense_SeqAIJ)); 4870 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4871 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJKron_C",MatSeqAIJKron_SeqAIJ)); 4872 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_SeqAIJ)); 4873 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetValuesCOO_C",MatSetValuesCOO_SeqAIJ)); 4874 PetscCall(MatCreate_SeqAIJ_Inode(B)); 4875 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4876 PetscCall(MatSeqAIJSetTypeFromOptions(B)); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4877 PetscFunctionReturn(0); 4878 } 4879 4880 /* 4881 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4882 */ 4883 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4884 { 4885 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data,*a = (Mat_SeqAIJ*)A->data; 4886 PetscInt m = A->rmap->n,i; 4887 4888 PetscFunctionBegin; 4889 PetscCheck(A->assembled || cpvalues == MAT_DO_NOT_COPY_VALUES,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot duplicate unassembled matrix"); 4890 4891 C->factortype = A->factortype; 4892 c->row = NULL; 4893 c->col = NULL; 4894 c->icol = NULL; 4895 c->reallocs = 0; 4896 4897 C->assembled = A->assembled; 4898 C->preallocated = A->preallocated; 4899 4900 if (A->preallocated) { 4901 PetscCall(PetscLayoutReference(A->rmap,&C->rmap)); 4902 PetscCall(PetscLayoutReference(A->cmap,&C->cmap)); 4903 4904 PetscCall(PetscMalloc1(m,&c->imax)); 4905 PetscCall(PetscMemcpy(c->imax,a->imax,m*sizeof(PetscInt))); 4906 PetscCall(PetscMalloc1(m,&c->ilen)); 4907 PetscCall(PetscMemcpy(c->ilen,a->ilen,m*sizeof(PetscInt))); 4908 PetscCall(PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt))); 4909 4910 /* allocate the matrix space */ 4911 if (mallocmatspace) { 4912 PetscCall(PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i)); 4913 PetscCall(PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt))); 4914 4915 c->singlemalloc = PETSC_TRUE; 4916 4917 PetscCall(PetscArraycpy(c->i,a->i,m+1)); 4918 if (m > 0) { 4919 PetscCall(PetscArraycpy(c->j,a->j,a->i[m])); 4920 if (cpvalues == MAT_COPY_VALUES) { 4921 const PetscScalar *aa; 4922 4923 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4924 PetscCall(PetscArraycpy(c->a,aa,a->i[m])); 4925 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4926 } else { 4927 PetscCall(PetscArrayzero(c->a,a->i[m])); 4928 } 4929 } 4930 } 4931 4932 c->ignorezeroentries = a->ignorezeroentries; 4933 c->roworiented = a->roworiented; 4934 c->nonew = a->nonew; 4935 if (a->diag) { 4936 PetscCall(PetscMalloc1(m+1,&c->diag)); 4937 PetscCall(PetscMemcpy(c->diag,a->diag,m*sizeof(PetscInt))); 4938 PetscCall(PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt))); 4939 } else c->diag = NULL; 4940 4941 c->solve_work = NULL; 4942 c->saved_values = NULL; 4943 c->idiag = NULL; 4944 c->ssor_work = NULL; 4945 c->keepnonzeropattern = a->keepnonzeropattern; 4946 c->free_a = PETSC_TRUE; 4947 c->free_ij = PETSC_TRUE; 4948 4949 c->rmax = a->rmax; 4950 c->nz = a->nz; 4951 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4952 4953 c->compressedrow.use = a->compressedrow.use; 4954 c->compressedrow.nrows = a->compressedrow.nrows; 4955 if (a->compressedrow.use) { 4956 i = a->compressedrow.nrows; 4957 PetscCall(PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex)); 4958 PetscCall(PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1)); 4959 PetscCall(PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i)); 4960 } else { 4961 c->compressedrow.use = PETSC_FALSE; 4962 c->compressedrow.i = NULL; 4963 c->compressedrow.rindex = NULL; 4964 } 4965 c->nonzerorowcnt = a->nonzerorowcnt; 4966 C->nonzerostate = A->nonzerostate; 4967 4968 PetscCall(MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C)); 4969 } 4970 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist)); 4971 PetscFunctionReturn(0); 4972 } 4973 4974 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4975 { 4976 PetscFunctionBegin; 4977 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 4978 PetscCall(MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n)); 4979 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4980 PetscCall(MatSetBlockSizesFromMats(*B,A,A)); 4981 } 4982 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 4983 PetscCall(MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE)); 4984 PetscFunctionReturn(0); 4985 } 4986 4987 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4988 { 4989 PetscBool isbinary, ishdf5; 4990 4991 PetscFunctionBegin; 4992 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 4993 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 4994 /* force binary viewer to load .info file if it has not yet done so */ 4995 PetscCall(PetscViewerSetUp(viewer)); 4996 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 4997 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5)); 4998 if (isbinary) { 4999 PetscCall(MatLoad_SeqAIJ_Binary(newMat,viewer)); 5000 } else if (ishdf5) { 5001 #if defined(PETSC_HAVE_HDF5) 5002 PetscCall(MatLoad_AIJ_HDF5(newMat,viewer)); 5003 #else 5004 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 5005 #endif 5006 } else { 5007 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name); 5008 } 5009 PetscFunctionReturn(0); 5010 } 5011 5012 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 5013 { 5014 Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data; 5015 PetscInt header[4],*rowlens,M,N,nz,sum,rows,cols,i; 5016 5017 PetscFunctionBegin; 5018 PetscCall(PetscViewerSetUp(viewer)); 5019 5020 /* read in matrix header */ 5021 PetscCall(PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT)); 5022 PetscCheck(header[0] == MAT_FILE_CLASSID,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 5023 M = header[1]; N = header[2]; nz = header[3]; 5024 PetscCheck(M >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%" PetscInt_FMT ") in file is negative",M); 5025 PetscCheck(N >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%" PetscInt_FMT ") in file is negative",N); 5026 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 5027 5028 /* set block sizes from the viewer's .info file */ 5029 PetscCall(MatLoad_Binary_BlockSizes(mat,viewer)); 5030 /* set local and global sizes if not set already */ 5031 if (mat->rmap->n < 0) mat->rmap->n = M; 5032 if (mat->cmap->n < 0) mat->cmap->n = N; 5033 if (mat->rmap->N < 0) mat->rmap->N = M; 5034 if (mat->cmap->N < 0) mat->cmap->N = N; 5035 PetscCall(PetscLayoutSetUp(mat->rmap)); 5036 PetscCall(PetscLayoutSetUp(mat->cmap)); 5037 5038 /* check if the matrix sizes are correct */ 5039 PetscCall(MatGetSize(mat,&rows,&cols)); 5040 PetscCheck(M == rows && N == cols,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")",M,N,rows,cols); 5041 5042 /* read in row lengths */ 5043 PetscCall(PetscMalloc1(M,&rowlens)); 5044 PetscCall(PetscViewerBinaryRead(viewer,rowlens,M,NULL,PETSC_INT)); 5045 /* check if sum(rowlens) is same as nz */ 5046 sum = 0; for (i=0; i<M; i++) sum += rowlens[i]; 5047 PetscCheck(sum == nz,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Inconsistent matrix data in file: nonzeros = %" PetscInt_FMT ", sum-row-lengths = %" PetscInt_FMT,nz,sum); 5048 /* preallocate and check sizes */ 5049 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat,0,rowlens)); 5050 PetscCall(MatGetSize(mat,&rows,&cols)); 5051 PetscCheck(M == rows && N == cols,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")",M,N,rows,cols); 5052 /* store row lengths */ 5053 PetscCall(PetscArraycpy(a->ilen,rowlens,M)); 5054 PetscCall(PetscFree(rowlens)); 5055 5056 /* fill in "i" row pointers */ 5057 a->i[0] = 0; for (i=0; i<M; i++) a->i[i+1] = a->i[i] + a->ilen[i]; 5058 /* read in "j" column indices */ 5059 PetscCall(PetscViewerBinaryRead(viewer,a->j,nz,NULL,PETSC_INT)); 5060 /* read in "a" nonzero values */ 5061 PetscCall(PetscViewerBinaryRead(viewer,a->a,nz,NULL,PETSC_SCALAR)); 5062 5063 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 5064 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 5065 PetscFunctionReturn(0); 5066 } 5067 5068 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 5069 { 5070 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 5071 const PetscScalar *aa,*ba; 5072 #if defined(PETSC_USE_COMPLEX) 5073 PetscInt k; 5074 #endif 5075 5076 PetscFunctionBegin; 5077 /* If the matrix dimensions are not equal,or no of nonzeros */ 5078 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 5079 *flg = PETSC_FALSE; 5080 PetscFunctionReturn(0); 5081 } 5082 5083 /* if the a->i are the same */ 5084 PetscCall(PetscArraycmp(a->i,b->i,A->rmap->n+1,flg)); 5085 if (!*flg) PetscFunctionReturn(0); 5086 5087 /* if a->j are the same */ 5088 PetscCall(PetscArraycmp(a->j,b->j,a->nz,flg)); 5089 if (!*flg) PetscFunctionReturn(0); 5090 5091 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 5092 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 5093 /* if a->a are the same */ 5094 #if defined(PETSC_USE_COMPLEX) 5095 for (k=0; k<a->nz; k++) { 5096 if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) { 5097 *flg = PETSC_FALSE; 5098 PetscFunctionReturn(0); 5099 } 5100 } 5101 #else 5102 PetscCall(PetscArraycmp(aa,ba,a->nz,flg)); 5103 #endif 5104 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 5105 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 5106 PetscFunctionReturn(0); 5107 } 5108 5109 /*@ 5110 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 5111 provided by the user. 5112 5113 Collective 5114 5115 Input Parameters: 5116 + comm - must be an MPI communicator of size 1 5117 . m - number of rows 5118 . n - number of columns 5119 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 5120 . j - column indices 5121 - a - matrix values 5122 5123 Output Parameter: 5124 . mat - the matrix 5125 5126 Level: intermediate 5127 5128 Notes: 5129 The i, j, and a arrays are not copied by this routine, the user must free these arrays 5130 once the matrix is destroyed and not before 5131 5132 You cannot set new nonzero locations into this matrix, that will generate an error. 5133 5134 The i and j indices are 0 based 5135 5136 The format which is used for the sparse matrix input, is equivalent to a 5137 row-major ordering.. i.e for the following matrix, the input data expected is 5138 as shown 5139 5140 $ 1 0 0 5141 $ 2 0 3 5142 $ 4 5 6 5143 $ 5144 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 5145 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 5146 $ v = {1,2,3,4,5,6} [size = 6] 5147 5148 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateMPIAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()` 5149 5150 @*/ 5151 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 5152 { 5153 PetscInt ii; 5154 Mat_SeqAIJ *aij; 5155 PetscInt jj; 5156 5157 PetscFunctionBegin; 5158 PetscCheck(m <= 0 || i[0] == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 5159 PetscCall(MatCreate(comm,mat)); 5160 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5161 /* PetscCall(MatSetBlockSizes(*mat,,)); */ 5162 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5163 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,NULL)); 5164 aij = (Mat_SeqAIJ*)(*mat)->data; 5165 PetscCall(PetscMalloc1(m,&aij->imax)); 5166 PetscCall(PetscMalloc1(m,&aij->ilen)); 5167 5168 aij->i = i; 5169 aij->j = j; 5170 aij->a = a; 5171 aij->singlemalloc = PETSC_FALSE; 5172 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 5173 aij->free_a = PETSC_FALSE; 5174 aij->free_ij = PETSC_FALSE; 5175 5176 for (ii=0,aij->nonzerorowcnt=0,aij->rmax=0; ii<m; ii++) { 5177 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 5178 if (PetscDefined(USE_DEBUG)) { 5179 PetscCheck(i[ii+1] - i[ii] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %" PetscInt_FMT " length = %" PetscInt_FMT,ii,i[ii+1] - i[ii]); 5180 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 5181 PetscCheck(j[jj] >= j[jj-1],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is not sorted",jj-i[ii],j[jj],ii); 5182 PetscCheck(j[jj] != j[jj-1],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is identical to previous entry",jj-i[ii],j[jj],ii); 5183 } 5184 } 5185 } 5186 if (PetscDefined(USE_DEBUG)) { 5187 for (ii=0; ii<aij->i[m]; ii++) { 5188 PetscCheck(j[ii] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 5189 PetscCheck(j[ii] <= n - 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 5190 } 5191 } 5192 5193 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5194 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5195 PetscFunctionReturn(0); 5196 } 5197 5198 /*@ 5199 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 5200 provided by the user. 5201 5202 Collective 5203 5204 Input Parameters: 5205 + comm - must be an MPI communicator of size 1 5206 . m - number of rows 5207 . n - number of columns 5208 . i - row indices 5209 . j - column indices 5210 . a - matrix values 5211 . nz - number of nonzeros 5212 - idx - if the i and j indices start with 1 use PETSC_TRUE otherwise use PETSC_FALSE 5213 5214 Output Parameter: 5215 . mat - the matrix 5216 5217 Level: intermediate 5218 5219 Example: 5220 For the following matrix, the input data expected is as shown (using 0 based indexing) 5221 .vb 5222 1 0 0 5223 2 0 3 5224 4 5 6 5225 5226 i = {0,1,1,2,2,2} 5227 j = {0,0,2,0,1,2} 5228 v = {1,2,3,4,5,6} 5229 .ve 5230 5231 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateSeqAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()`, `MatSetValuesCOO()` 5232 5233 @*/ 5234 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 5235 { 5236 PetscInt ii, *nnz, one = 1,row,col; 5237 5238 PetscFunctionBegin; 5239 PetscCall(PetscCalloc1(m,&nnz)); 5240 for (ii = 0; ii < nz; ii++) { 5241 nnz[i[ii] - !!idx] += 1; 5242 } 5243 PetscCall(MatCreate(comm,mat)); 5244 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5245 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5246 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz)); 5247 for (ii = 0; ii < nz; ii++) { 5248 if (idx) { 5249 row = i[ii] - 1; 5250 col = j[ii] - 1; 5251 } else { 5252 row = i[ii]; 5253 col = j[ii]; 5254 } 5255 PetscCall(MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES)); 5256 } 5257 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5258 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5259 PetscCall(PetscFree(nnz)); 5260 PetscFunctionReturn(0); 5261 } 5262 5263 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 5264 { 5265 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 5266 5267 PetscFunctionBegin; 5268 a->idiagvalid = PETSC_FALSE; 5269 a->ibdiagvalid = PETSC_FALSE; 5270 5271 PetscCall(MatSeqAIJInvalidateDiagonal_Inode(A)); 5272 PetscFunctionReturn(0); 5273 } 5274 5275 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 5276 { 5277 PetscMPIInt size; 5278 5279 PetscFunctionBegin; 5280 PetscCallMPI(MPI_Comm_size(comm,&size)); 5281 if (size == 1) { 5282 if (scall == MAT_INITIAL_MATRIX) { 5283 PetscCall(MatDuplicate(inmat,MAT_COPY_VALUES,outmat)); 5284 } else { 5285 PetscCall(MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN)); 5286 } 5287 } else { 5288 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat)); 5289 } 5290 PetscFunctionReturn(0); 5291 } 5292 5293 /* 5294 Permute A into C's *local* index space using rowemb,colemb. 5295 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 5296 of [0,m), colemb is in [0,n). 5297 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 5298 */ 5299 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 5300 { 5301 /* If making this function public, change the error returned in this function away from _PLIB. */ 5302 Mat_SeqAIJ *Baij; 5303 PetscBool seqaij; 5304 PetscInt m,n,*nz,i,j,count; 5305 PetscScalar v; 5306 const PetscInt *rowindices,*colindices; 5307 5308 PetscFunctionBegin; 5309 if (!B) PetscFunctionReturn(0); 5310 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 5311 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij)); 5312 PetscCheck(seqaij,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 5313 if (rowemb) { 5314 PetscCall(ISGetLocalSize(rowemb,&m)); 5315 PetscCheck(m == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Row IS of size %" PetscInt_FMT " is incompatible with matrix row size %" PetscInt_FMT,m,B->rmap->n); 5316 } else { 5317 PetscCheck(C->rmap->n == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 5318 } 5319 if (colemb) { 5320 PetscCall(ISGetLocalSize(colemb,&n)); 5321 PetscCheck(n == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diag col IS of size %" PetscInt_FMT " is incompatible with input matrix col size %" PetscInt_FMT,n,B->cmap->n); 5322 } else { 5323 PetscCheck(C->cmap->n == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 5324 } 5325 5326 Baij = (Mat_SeqAIJ*)(B->data); 5327 if (pattern == DIFFERENT_NONZERO_PATTERN) { 5328 PetscCall(PetscMalloc1(B->rmap->n,&nz)); 5329 for (i=0; i<B->rmap->n; i++) { 5330 nz[i] = Baij->i[i+1] - Baij->i[i]; 5331 } 5332 PetscCall(MatSeqAIJSetPreallocation(C,0,nz)); 5333 PetscCall(PetscFree(nz)); 5334 } 5335 if (pattern == SUBSET_NONZERO_PATTERN) { 5336 PetscCall(MatZeroEntries(C)); 5337 } 5338 count = 0; 5339 rowindices = NULL; 5340 colindices = NULL; 5341 if (rowemb) { 5342 PetscCall(ISGetIndices(rowemb,&rowindices)); 5343 } 5344 if (colemb) { 5345 PetscCall(ISGetIndices(colemb,&colindices)); 5346 } 5347 for (i=0; i<B->rmap->n; i++) { 5348 PetscInt row; 5349 row = i; 5350 if (rowindices) row = rowindices[i]; 5351 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 5352 PetscInt col; 5353 col = Baij->j[count]; 5354 if (colindices) col = colindices[col]; 5355 v = Baij->a[count]; 5356 PetscCall(MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES)); 5357 ++count; 5358 } 5359 } 5360 /* FIXME: set C's nonzerostate correctly. */ 5361 /* Assembly for C is necessary. */ 5362 C->preallocated = PETSC_TRUE; 5363 C->assembled = PETSC_TRUE; 5364 C->was_assembled = PETSC_FALSE; 5365 PetscFunctionReturn(0); 5366 } 5367 5368 PetscFunctionList MatSeqAIJList = NULL; 5369 5370 /*@C 5371 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 5372 5373 Collective on Mat 5374 5375 Input Parameters: 5376 + mat - the matrix object 5377 - matype - matrix type 5378 5379 Options Database Key: 5380 . -mat_seqai_type <method> - for example seqaijcrl 5381 5382 Level: intermediate 5383 5384 .seealso: `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`, `Mat` 5385 @*/ 5386 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 5387 { 5388 PetscBool sametype; 5389 PetscErrorCode (*r)(Mat,MatType,MatReuse,Mat*); 5390 5391 PetscFunctionBegin; 5392 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 5393 PetscCall(PetscObjectTypeCompare((PetscObject)mat,matype,&sametype)); 5394 if (sametype) PetscFunctionReturn(0); 5395 5396 PetscCall(PetscFunctionListFind(MatSeqAIJList,matype,&r)); 5397 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 5398 PetscCall((*r)(mat,matype,MAT_INPLACE_MATRIX,&mat)); 5399 PetscFunctionReturn(0); 5400 } 5401 5402 /*@C 5403 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 5404 5405 Not Collective 5406 5407 Input Parameters: 5408 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 5409 - function - routine to convert to subtype 5410 5411 Notes: 5412 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 5413 5414 Then, your matrix can be chosen with the procedural interface at runtime via the option 5415 $ -mat_seqaij_type my_mat 5416 5417 Level: advanced 5418 5419 .seealso: `MatSeqAIJRegisterAll()` 5420 5421 Level: advanced 5422 @*/ 5423 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 5424 { 5425 PetscFunctionBegin; 5426 PetscCall(MatInitializePackage()); 5427 PetscCall(PetscFunctionListAdd(&MatSeqAIJList,sname,function)); 5428 PetscFunctionReturn(0); 5429 } 5430 5431 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 5432 5433 /*@C 5434 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 5435 5436 Not Collective 5437 5438 Level: advanced 5439 5440 .seealso: `MatRegisterAll()`, `MatSeqAIJRegister()` 5441 @*/ 5442 PetscErrorCode MatSeqAIJRegisterAll(void) 5443 { 5444 PetscFunctionBegin; 5445 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 5446 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 5447 5448 PetscCall(MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL)); 5449 PetscCall(MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM)); 5450 PetscCall(MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL)); 5451 #if defined(PETSC_HAVE_MKL_SPARSE) 5452 PetscCall(MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL)); 5453 #endif 5454 #if defined(PETSC_HAVE_CUDA) 5455 PetscCall(MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 5456 #endif 5457 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 5458 PetscCall(MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos)); 5459 #endif 5460 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 5461 PetscCall(MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL)); 5462 #endif 5463 PetscFunctionReturn(0); 5464 } 5465 5466 /* 5467 Special version for direct calls from Fortran 5468 */ 5469 #include <petsc/private/fortranimpl.h> 5470 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5471 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 5472 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 5473 #define matsetvaluesseqaij_ matsetvaluesseqaij 5474 #endif 5475 5476 /* Change these macros so can be used in void function */ 5477 5478 /* Change these macros so can be used in void function */ 5479 /* Identical to PetscCallVoid, except it assigns to *_ierr */ 5480 #undef PetscCall 5481 #define PetscCall(...) do { \ 5482 PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \ 5483 if (PetscUnlikely(ierr_msv_mpiaij)) { \ 5484 *_ierr = PetscError(PETSC_COMM_SELF,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr_msv_mpiaij,PETSC_ERROR_REPEAT," "); \ 5485 return; \ 5486 } \ 5487 } while (0) 5488 5489 #undef SETERRQ 5490 #define SETERRQ(comm,ierr,...) do { \ 5491 *_ierr = PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr,PETSC_ERROR_INITIAL,__VA_ARGS__); \ 5492 return; \ 5493 } while (0) 5494 5495 PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 5496 { 5497 Mat A = *AA; 5498 PetscInt m = *mm, n = *nn; 5499 InsertMode is = *isis; 5500 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 5501 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 5502 PetscInt *imax,*ai,*ailen; 5503 PetscInt *aj,nonew = a->nonew,lastcol = -1; 5504 MatScalar *ap,value,*aa; 5505 PetscBool ignorezeroentries = a->ignorezeroentries; 5506 PetscBool roworiented = a->roworiented; 5507 5508 PetscFunctionBegin; 5509 MatCheckPreallocated(A,1); 5510 imax = a->imax; 5511 ai = a->i; 5512 ailen = a->ilen; 5513 aj = a->j; 5514 aa = a->a; 5515 5516 for (k=0; k<m; k++) { /* loop over added rows */ 5517 row = im[k]; 5518 if (row < 0) continue; 5519 PetscCheck(row < A->rmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 5520 rp = aj + ai[row]; ap = aa + ai[row]; 5521 rmax = imax[row]; nrow = ailen[row]; 5522 low = 0; 5523 high = nrow; 5524 for (l=0; l<n; l++) { /* loop over added columns */ 5525 if (in[l] < 0) continue; 5526 PetscCheck(in[l] < A->cmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 5527 col = in[l]; 5528 if (roworiented) value = v[l + k*n]; 5529 else value = v[k + l*m]; 5530 5531 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 5532 5533 if (col <= lastcol) low = 0; 5534 else high = nrow; 5535 lastcol = col; 5536 while (high-low > 5) { 5537 t = (low+high)/2; 5538 if (rp[t] > col) high = t; 5539 else low = t; 5540 } 5541 for (i=low; i<high; i++) { 5542 if (rp[i] > col) break; 5543 if (rp[i] == col) { 5544 if (is == ADD_VALUES) ap[i] += value; 5545 else ap[i] = value; 5546 goto noinsert; 5547 } 5548 } 5549 if (value == 0.0 && ignorezeroentries) goto noinsert; 5550 if (nonew == 1) goto noinsert; 5551 PetscCheck(nonew != -1,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 5552 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 5553 N = nrow++ - 1; a->nz++; high++; 5554 /* shift up all the later entries in this row */ 5555 for (ii=N; ii>=i; ii--) { 5556 rp[ii+1] = rp[ii]; 5557 ap[ii+1] = ap[ii]; 5558 } 5559 rp[i] = col; 5560 ap[i] = value; 5561 A->nonzerostate++; 5562 noinsert:; 5563 low = i + 1; 5564 } 5565 ailen[row] = nrow; 5566 } 5567 PetscFunctionReturnVoid(); 5568 } 5569 /* Undefining these here since they were redefined from their original definition above! No 5570 * other PETSc functions should be defined past this point, as it is impossible to recover the 5571 * original definitions */ 5572 #undef PetscCall 5573 #undef SETERRQ 5574