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 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 MatCreateGraph_Simple_AIJ, 3633 MatFilter_AIJ 3634 }; 3635 3636 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3637 { 3638 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3639 PetscInt i,nz,n; 3640 3641 PetscFunctionBegin; 3642 nz = aij->maxnz; 3643 n = mat->rmap->n; 3644 for (i=0; i<nz; i++) { 3645 aij->j[i] = indices[i]; 3646 } 3647 aij->nz = nz; 3648 for (i=0; i<n; i++) { 3649 aij->ilen[i] = aij->imax[i]; 3650 } 3651 PetscFunctionReturn(0); 3652 } 3653 3654 /* 3655 * Given a sparse matrix with global column indices, compact it by using a local column space. 3656 * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3657 */ 3658 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3659 { 3660 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3661 PetscTable gid1_lid1; 3662 PetscTablePosition tpos; 3663 PetscInt gid,lid,i,ec,nz = aij->nz; 3664 PetscInt *garray,*jj = aij->j; 3665 3666 PetscFunctionBegin; 3667 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3668 PetscValidPointer(mapping,2); 3669 /* use a table */ 3670 PetscCall(PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1)); 3671 ec = 0; 3672 for (i=0; i<nz; i++) { 3673 PetscInt data,gid1 = jj[i] + 1; 3674 PetscCall(PetscTableFind(gid1_lid1,gid1,&data)); 3675 if (!data) { 3676 /* one based table */ 3677 PetscCall(PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES)); 3678 } 3679 } 3680 /* form array of columns we need */ 3681 PetscCall(PetscMalloc1(ec,&garray)); 3682 PetscCall(PetscTableGetHeadPosition(gid1_lid1,&tpos)); 3683 while (tpos) { 3684 PetscCall(PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid)); 3685 gid--; 3686 lid--; 3687 garray[lid] = gid; 3688 } 3689 PetscCall(PetscSortInt(ec,garray)); /* sort, and rebuild */ 3690 PetscCall(PetscTableRemoveAll(gid1_lid1)); 3691 for (i=0; i<ec; i++) { 3692 PetscCall(PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES)); 3693 } 3694 /* compact out the extra columns in B */ 3695 for (i=0; i<nz; i++) { 3696 PetscInt gid1 = jj[i] + 1; 3697 PetscCall(PetscTableFind(gid1_lid1,gid1,&lid)); 3698 lid--; 3699 jj[i] = lid; 3700 } 3701 PetscCall(PetscLayoutDestroy(&mat->cmap)); 3702 PetscCall(PetscTableDestroy(&gid1_lid1)); 3703 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat),ec,ec,1,&mat->cmap)); 3704 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping)); 3705 PetscCall(ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH)); 3706 PetscFunctionReturn(0); 3707 } 3708 3709 /*@ 3710 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3711 in the matrix. 3712 3713 Input Parameters: 3714 + mat - the SeqAIJ matrix 3715 - indices - the column indices 3716 3717 Level: advanced 3718 3719 Notes: 3720 This can be called if you have precomputed the nonzero structure of the 3721 matrix and want to provide it to the matrix object to improve the performance 3722 of the MatSetValues() operation. 3723 3724 You MUST have set the correct numbers of nonzeros per row in the call to 3725 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3726 3727 MUST be called before any calls to MatSetValues(); 3728 3729 The indices should start with zero, not one. 3730 3731 @*/ 3732 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3733 { 3734 PetscFunctionBegin; 3735 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3736 PetscValidIntPointer(indices,2); 3737 PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices)); 3738 PetscFunctionReturn(0); 3739 } 3740 3741 /* ----------------------------------------------------------------------------------------*/ 3742 3743 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3744 { 3745 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3746 size_t nz = aij->i[mat->rmap->n]; 3747 3748 PetscFunctionBegin; 3749 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3750 3751 /* allocate space for values if not already there */ 3752 if (!aij->saved_values) { 3753 PetscCall(PetscMalloc1(nz+1,&aij->saved_values)); 3754 PetscCall(PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar))); 3755 } 3756 3757 /* copy values over */ 3758 PetscCall(PetscArraycpy(aij->saved_values,aij->a,nz)); 3759 PetscFunctionReturn(0); 3760 } 3761 3762 /*@ 3763 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3764 example, reuse of the linear part of a Jacobian, while recomputing the 3765 nonlinear portion. 3766 3767 Collect on Mat 3768 3769 Input Parameters: 3770 . mat - the matrix (currently only AIJ matrices support this option) 3771 3772 Level: advanced 3773 3774 Common Usage, with SNESSolve(): 3775 $ Create Jacobian matrix 3776 $ Set linear terms into matrix 3777 $ Apply boundary conditions to matrix, at this time matrix must have 3778 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3779 $ boundary conditions again will not change the nonzero structure 3780 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3781 $ ierr = MatStoreValues(mat); 3782 $ Call SNESSetJacobian() with matrix 3783 $ In your Jacobian routine 3784 $ ierr = MatRetrieveValues(mat); 3785 $ Set nonlinear terms in matrix 3786 3787 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3788 $ // build linear portion of Jacobian 3789 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3790 $ ierr = MatStoreValues(mat); 3791 $ loop over nonlinear iterations 3792 $ ierr = MatRetrieveValues(mat); 3793 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3794 $ // call MatAssemblyBegin/End() on matrix 3795 $ Solve linear system with Jacobian 3796 $ endloop 3797 3798 Notes: 3799 Matrix must already be assemblied before calling this routine 3800 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3801 calling this routine. 3802 3803 When this is called multiple times it overwrites the previous set of stored values 3804 and does not allocated additional space. 3805 3806 .seealso: `MatRetrieveValues()` 3807 3808 @*/ 3809 PetscErrorCode MatStoreValues(Mat mat) 3810 { 3811 PetscFunctionBegin; 3812 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3813 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3814 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3815 PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat)); 3816 PetscFunctionReturn(0); 3817 } 3818 3819 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3820 { 3821 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3822 PetscInt nz = aij->i[mat->rmap->n]; 3823 3824 PetscFunctionBegin; 3825 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3826 PetscCheck(aij->saved_values,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3827 /* copy values over */ 3828 PetscCall(PetscArraycpy(aij->a,aij->saved_values,nz)); 3829 PetscFunctionReturn(0); 3830 } 3831 3832 /*@ 3833 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3834 example, reuse of the linear part of a Jacobian, while recomputing the 3835 nonlinear portion. 3836 3837 Collect on Mat 3838 3839 Input Parameters: 3840 . mat - the matrix (currently only AIJ matrices support this option) 3841 3842 Level: advanced 3843 3844 .seealso: `MatStoreValues()` 3845 3846 @*/ 3847 PetscErrorCode MatRetrieveValues(Mat mat) 3848 { 3849 PetscFunctionBegin; 3850 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3851 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3852 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3853 PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat)); 3854 PetscFunctionReturn(0); 3855 } 3856 3857 /* --------------------------------------------------------------------------------*/ 3858 /*@C 3859 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3860 (the default parallel PETSc format). For good matrix assembly performance 3861 the user should preallocate the matrix storage by setting the parameter nz 3862 (or the array nnz). By setting these parameters accurately, performance 3863 during matrix assembly can be increased by more than a factor of 50. 3864 3865 Collective 3866 3867 Input Parameters: 3868 + comm - MPI communicator, set to PETSC_COMM_SELF 3869 . m - number of rows 3870 . n - number of columns 3871 . nz - number of nonzeros per row (same for all rows) 3872 - nnz - array containing the number of nonzeros in the various rows 3873 (possibly different for each row) or NULL 3874 3875 Output Parameter: 3876 . A - the matrix 3877 3878 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3879 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3880 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3881 3882 Notes: 3883 If nnz is given then nz is ignored 3884 3885 The AIJ format (also called the Yale sparse matrix format or 3886 compressed row storage), is fully compatible with standard Fortran 77 3887 storage. That is, the stored row and column indices can begin at 3888 either one (as in Fortran) or zero. See the users' manual for details. 3889 3890 Specify the preallocated storage with either nz or nnz (not both). 3891 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3892 allocation. For large problems you MUST preallocate memory or you 3893 will get TERRIBLE performance, see the users' manual chapter on matrices. 3894 3895 By default, this format uses inodes (identical nodes) when possible, to 3896 improve numerical efficiency of matrix-vector products and solves. We 3897 search for consecutive rows with the same nonzero structure, thereby 3898 reusing matrix information to achieve increased efficiency. 3899 3900 Options Database Keys: 3901 + -mat_no_inode - Do not use inodes 3902 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3903 3904 Level: intermediate 3905 3906 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()` 3907 3908 @*/ 3909 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3910 { 3911 PetscFunctionBegin; 3912 PetscCall(MatCreate(comm,A)); 3913 PetscCall(MatSetSizes(*A,m,n,m,n)); 3914 PetscCall(MatSetType(*A,MATSEQAIJ)); 3915 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz)); 3916 PetscFunctionReturn(0); 3917 } 3918 3919 /*@C 3920 MatSeqAIJSetPreallocation - For good matrix assembly performance 3921 the user should preallocate the matrix storage by setting the parameter nz 3922 (or the array nnz). By setting these parameters accurately, performance 3923 during matrix assembly can be increased by more than a factor of 50. 3924 3925 Collective 3926 3927 Input Parameters: 3928 + B - The matrix 3929 . nz - number of nonzeros per row (same for all rows) 3930 - nnz - array containing the number of nonzeros in the various rows 3931 (possibly different for each row) or NULL 3932 3933 Notes: 3934 If nnz is given then nz is ignored 3935 3936 The AIJ format (also called the Yale sparse matrix format or 3937 compressed row storage), is fully compatible with standard Fortran 77 3938 storage. That is, the stored row and column indices can begin at 3939 either one (as in Fortran) or zero. See the users' manual for details. 3940 3941 Specify the preallocated storage with either nz or nnz (not both). 3942 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3943 allocation. For large problems you MUST preallocate memory or you 3944 will get TERRIBLE performance, see the users' manual chapter on matrices. 3945 3946 You can call MatGetInfo() to get information on how effective the preallocation was; 3947 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3948 You can also run with the option -info and look for messages with the string 3949 malloc in them to see if additional memory allocation was needed. 3950 3951 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3952 entries or columns indices 3953 3954 By default, this format uses inodes (identical nodes) when possible, to 3955 improve numerical efficiency of matrix-vector products and solves. We 3956 search for consecutive rows with the same nonzero structure, thereby 3957 reusing matrix information to achieve increased efficiency. 3958 3959 Options Database Keys: 3960 + -mat_no_inode - Do not use inodes 3961 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3962 3963 Level: intermediate 3964 3965 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()`, `MatGetInfo()`, 3966 `MatSeqAIJSetTotalPreallocation()` 3967 3968 @*/ 3969 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3970 { 3971 PetscFunctionBegin; 3972 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3973 PetscValidType(B,1); 3974 PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz)); 3975 PetscFunctionReturn(0); 3976 } 3977 3978 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3979 { 3980 Mat_SeqAIJ *b; 3981 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3982 PetscInt i; 3983 3984 PetscFunctionBegin; 3985 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3986 if (nz == MAT_SKIP_ALLOCATION) { 3987 skipallocation = PETSC_TRUE; 3988 nz = 0; 3989 } 3990 PetscCall(PetscLayoutSetUp(B->rmap)); 3991 PetscCall(PetscLayoutSetUp(B->cmap)); 3992 3993 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3994 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %" PetscInt_FMT,nz); 3995 if (PetscUnlikelyDebug(nnz)) { 3996 for (i=0; i<B->rmap->n; i++) { 3997 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]); 3998 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); 3999 } 4000 } 4001 4002 B->preallocated = PETSC_TRUE; 4003 4004 b = (Mat_SeqAIJ*)B->data; 4005 4006 if (!skipallocation) { 4007 if (!b->imax) { 4008 PetscCall(PetscMalloc1(B->rmap->n,&b->imax)); 4009 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4010 } 4011 if (!b->ilen) { 4012 /* b->ilen will count nonzeros in each row so far. */ 4013 PetscCall(PetscCalloc1(B->rmap->n,&b->ilen)); 4014 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4015 } else { 4016 PetscCall(PetscMemzero(b->ilen,B->rmap->n*sizeof(PetscInt))); 4017 } 4018 if (!b->ipre) { 4019 PetscCall(PetscMalloc1(B->rmap->n,&b->ipre)); 4020 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4021 } 4022 if (!nnz) { 4023 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 4024 else if (nz < 0) nz = 1; 4025 nz = PetscMin(nz,B->cmap->n); 4026 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 4027 nz = nz*B->rmap->n; 4028 } else { 4029 PetscInt64 nz64 = 0; 4030 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 4031 PetscCall(PetscIntCast(nz64,&nz)); 4032 } 4033 4034 /* allocate the matrix space */ 4035 /* FIXME: should B's old memory be unlogged? */ 4036 PetscCall(MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i)); 4037 if (B->structure_only) { 4038 PetscCall(PetscMalloc1(nz,&b->j)); 4039 PetscCall(PetscMalloc1(B->rmap->n+1,&b->i)); 4040 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt))); 4041 } else { 4042 PetscCall(PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i)); 4043 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)))); 4044 } 4045 b->i[0] = 0; 4046 for (i=1; i<B->rmap->n+1; i++) { 4047 b->i[i] = b->i[i-1] + b->imax[i-1]; 4048 } 4049 if (B->structure_only) { 4050 b->singlemalloc = PETSC_FALSE; 4051 b->free_a = PETSC_FALSE; 4052 } else { 4053 b->singlemalloc = PETSC_TRUE; 4054 b->free_a = PETSC_TRUE; 4055 } 4056 b->free_ij = PETSC_TRUE; 4057 } else { 4058 b->free_a = PETSC_FALSE; 4059 b->free_ij = PETSC_FALSE; 4060 } 4061 4062 if (b->ipre && nnz != b->ipre && b->imax) { 4063 /* reserve user-requested sparsity */ 4064 PetscCall(PetscArraycpy(b->ipre,b->imax,B->rmap->n)); 4065 } 4066 4067 b->nz = 0; 4068 b->maxnz = nz; 4069 B->info.nz_unneeded = (double)b->maxnz; 4070 if (realalloc) PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE)); 4071 B->was_assembled = PETSC_FALSE; 4072 B->assembled = PETSC_FALSE; 4073 /* We simply deem preallocation has changed nonzero state. Updating the state 4074 will give clients (like AIJKokkos) a chance to know something has happened. 4075 */ 4076 B->nonzerostate++; 4077 PetscFunctionReturn(0); 4078 } 4079 4080 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 4081 { 4082 Mat_SeqAIJ *a; 4083 PetscInt i; 4084 4085 PetscFunctionBegin; 4086 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4087 4088 /* Check local size. If zero, then return */ 4089 if (!A->rmap->n) PetscFunctionReturn(0); 4090 4091 a = (Mat_SeqAIJ*)A->data; 4092 /* if no saved info, we error out */ 4093 PetscCheck(a->ipre,PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"No saved preallocation info "); 4094 4095 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 "); 4096 4097 PetscCall(PetscArraycpy(a->imax,a->ipre,A->rmap->n)); 4098 PetscCall(PetscArrayzero(a->ilen,A->rmap->n)); 4099 a->i[0] = 0; 4100 for (i=1; i<A->rmap->n+1; i++) { 4101 a->i[i] = a->i[i-1] + a->imax[i-1]; 4102 } 4103 A->preallocated = PETSC_TRUE; 4104 a->nz = 0; 4105 a->maxnz = a->i[A->rmap->n]; 4106 A->info.nz_unneeded = (double)a->maxnz; 4107 A->was_assembled = PETSC_FALSE; 4108 A->assembled = PETSC_FALSE; 4109 PetscFunctionReturn(0); 4110 } 4111 4112 /*@ 4113 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 4114 4115 Input Parameters: 4116 + B - the matrix 4117 . i - the indices into j for the start of each row (starts with zero) 4118 . j - the column indices for each row (starts with zero) these must be sorted for each row 4119 - v - optional values in the matrix 4120 4121 Level: developer 4122 4123 Notes: 4124 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 4125 4126 This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero 4127 structure will be the union of all the previous nonzero structures. 4128 4129 Developer Notes: 4130 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 4131 then just copies the v values directly with PetscMemcpy(). 4132 4133 This routine could also take a PetscCopyMode argument to allow sharing the values instead of always copying them. 4134 4135 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatSeqAIJSetPreallocation()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MatResetPreallocation()` 4136 @*/ 4137 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 4138 { 4139 PetscFunctionBegin; 4140 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 4141 PetscValidType(B,1); 4142 PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v)); 4143 PetscFunctionReturn(0); 4144 } 4145 4146 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 4147 { 4148 PetscInt i; 4149 PetscInt m,n; 4150 PetscInt nz; 4151 PetscInt *nnz; 4152 4153 PetscFunctionBegin; 4154 PetscCheck(Ii[0] == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %" PetscInt_FMT, Ii[0]); 4155 4156 PetscCall(PetscLayoutSetUp(B->rmap)); 4157 PetscCall(PetscLayoutSetUp(B->cmap)); 4158 4159 PetscCall(MatGetSize(B, &m, &n)); 4160 PetscCall(PetscMalloc1(m+1, &nnz)); 4161 for (i = 0; i < m; i++) { 4162 nz = Ii[i+1]- Ii[i]; 4163 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 4164 nnz[i] = nz; 4165 } 4166 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 4167 PetscCall(PetscFree(nnz)); 4168 4169 for (i = 0; i < m; i++) { 4170 PetscCall(MatSetValues_SeqAIJ(B, 1, &i, Ii[i+1] - Ii[i], J+Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES)); 4171 } 4172 4173 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 4174 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 4175 4176 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 4177 PetscFunctionReturn(0); 4178 } 4179 4180 /*@ 4181 MatSeqAIJKron - Computes C, the Kronecker product of A and B. 4182 4183 Input Parameters: 4184 + A - left-hand side matrix 4185 . B - right-hand side matrix 4186 - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 4187 4188 Output Parameter: 4189 . C - Kronecker product of A and B 4190 4191 Level: intermediate 4192 4193 Notes: 4194 MAT_REUSE_MATRIX can only be used when the nonzero structure of the product matrix has not changed from that last call to MatSeqAIJKron(). 4195 4196 .seealso: `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATKAIJ`, `MatReuse` 4197 @*/ 4198 PetscErrorCode MatSeqAIJKron(Mat A,Mat B,MatReuse reuse,Mat *C) 4199 { 4200 PetscFunctionBegin; 4201 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4202 PetscValidType(A,1); 4203 PetscValidHeaderSpecific(B,MAT_CLASSID,2); 4204 PetscValidType(B,2); 4205 PetscValidPointer(C,4); 4206 if (reuse == MAT_REUSE_MATRIX) { 4207 PetscValidHeaderSpecific(*C,MAT_CLASSID,4); 4208 PetscValidType(*C,4); 4209 } 4210 PetscTryMethod(A,"MatSeqAIJKron_C",(Mat,Mat,MatReuse,Mat*),(A,B,reuse,C)); 4211 PetscFunctionReturn(0); 4212 } 4213 4214 PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A,Mat B,MatReuse reuse,Mat *C) 4215 { 4216 Mat newmat; 4217 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4218 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 4219 PetscScalar *v; 4220 const PetscScalar *aa,*ba; 4221 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; 4222 PetscBool flg; 4223 4224 PetscFunctionBegin; 4225 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4226 PetscCheck(A->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4227 PetscCheck(!B->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4228 PetscCheck(B->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4229 PetscCall(PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&flg)); 4230 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatType %s",((PetscObject)B)->type_name); 4231 PetscCheck(reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatReuse %d",(int)reuse); 4232 if (reuse == MAT_INITIAL_MATRIX) { 4233 PetscCall(PetscMalloc2(am*bm+1,&i,a->i[am]*b->i[bm],&j)); 4234 PetscCall(MatCreate(PETSC_COMM_SELF,&newmat)); 4235 PetscCall(MatSetSizes(newmat,am*bm,an*bn,am*bm,an*bn)); 4236 PetscCall(MatSetType(newmat,MATAIJ)); 4237 i[0] = 0; 4238 for (m = 0; m < am; ++m) { 4239 for (p = 0; p < bm; ++p) { 4240 i[m*bm + p + 1] = i[m*bm + p] + (a->i[m+1] - a->i[m]) * (b->i[p+1] - b->i[p]); 4241 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4242 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4243 j[nnz++] = a->j[n]*bn + b->j[q]; 4244 } 4245 } 4246 } 4247 } 4248 PetscCall(MatSeqAIJSetPreallocationCSR(newmat,i,j,NULL)); 4249 *C = newmat; 4250 PetscCall(PetscFree2(i,j)); 4251 nnz = 0; 4252 } 4253 PetscCall(MatSeqAIJGetArray(*C,&v)); 4254 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4255 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 4256 for (m = 0; m < am; ++m) { 4257 for (p = 0; p < bm; ++p) { 4258 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4259 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4260 v[nnz++] = aa[n] * ba[q]; 4261 } 4262 } 4263 } 4264 } 4265 PetscCall(MatSeqAIJRestoreArray(*C,&v)); 4266 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 4267 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 4268 PetscFunctionReturn(0); 4269 } 4270 4271 #include <../src/mat/impls/dense/seq/dense.h> 4272 #include <petsc/private/kernels/petscaxpy.h> 4273 4274 /* 4275 Computes (B'*A')' since computing B*A directly is untenable 4276 4277 n p p 4278 [ ] [ ] [ ] 4279 m [ A ] * n [ B ] = m [ C ] 4280 [ ] [ ] [ ] 4281 4282 */ 4283 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 4284 { 4285 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 4286 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 4287 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 4288 PetscInt i,j,n,m,q,p; 4289 const PetscInt *ii,*idx; 4290 const PetscScalar *b,*a,*a_q; 4291 PetscScalar *c,*c_q; 4292 PetscInt clda = sub_c->lda; 4293 PetscInt alda = sub_a->lda; 4294 4295 PetscFunctionBegin; 4296 m = A->rmap->n; 4297 n = A->cmap->n; 4298 p = B->cmap->n; 4299 a = sub_a->v; 4300 b = sub_b->a; 4301 c = sub_c->v; 4302 if (clda == m) { 4303 PetscCall(PetscArrayzero(c,m*p)); 4304 } else { 4305 for (j=0;j<p;j++) 4306 for (i=0;i<m;i++) 4307 c[j*clda + i] = 0.0; 4308 } 4309 ii = sub_b->i; 4310 idx = sub_b->j; 4311 for (i=0; i<n; i++) { 4312 q = ii[i+1] - ii[i]; 4313 while (q-->0) { 4314 c_q = c + clda*(*idx); 4315 a_q = a + alda*i; 4316 PetscKernelAXPY(c_q,*b,a_q,m); 4317 idx++; 4318 b++; 4319 } 4320 } 4321 PetscFunctionReturn(0); 4322 } 4323 4324 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C) 4325 { 4326 PetscInt m=A->rmap->n,n=B->cmap->n; 4327 PetscBool cisdense; 4328 4329 PetscFunctionBegin; 4330 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); 4331 PetscCall(MatSetSizes(C,m,n,m,n)); 4332 PetscCall(MatSetBlockSizesFromMats(C,A,B)); 4333 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 4334 if (!cisdense) { 4335 PetscCall(MatSetType(C,MATDENSE)); 4336 } 4337 PetscCall(MatSetUp(C)); 4338 4339 C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4340 PetscFunctionReturn(0); 4341 } 4342 4343 /* ----------------------------------------------------------------*/ 4344 /*MC 4345 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4346 based on compressed sparse row format. 4347 4348 Options Database Keys: 4349 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4350 4351 Level: beginner 4352 4353 Notes: 4354 MatSetValues() may be called for this matrix type with a NULL argument for the numerical values, 4355 in this case the values associated with the rows and columns one passes in are set to zero 4356 in the matrix 4357 4358 MatSetOptions(,MAT_STRUCTURE_ONLY,PETSC_TRUE) may be called for this matrix type. In this no 4359 space is allocated for the nonzero entries and any entries passed with MatSetValues() are ignored 4360 4361 Developer Notes: 4362 It would be nice if all matrix formats supported passing NULL in for the numerical values 4363 4364 .seealso: `MatCreateSeqAIJ()`, `MatSetFromOptions()`, `MatSetType()`, `MatCreate()`, `MatType`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4365 M*/ 4366 4367 /*MC 4368 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4369 4370 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 4371 and MATMPIAIJ otherwise. As a result, for single process communicators, 4372 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4373 for communicators controlling multiple processes. It is recommended that you call both of 4374 the above preallocation routines for simplicity. 4375 4376 Options Database Keys: 4377 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 4378 4379 Developer Notes: 4380 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 4381 enough exist. 4382 4383 Level: beginner 4384 4385 .seealso: `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATMPIAIJ`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4386 M*/ 4387 4388 /*MC 4389 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4390 4391 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 4392 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 4393 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4394 for communicators controlling multiple processes. It is recommended that you call both of 4395 the above preallocation routines for simplicity. 4396 4397 Options Database Keys: 4398 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 4399 4400 Level: beginner 4401 4402 .seealso: `MatCreateMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL` 4403 M*/ 4404 4405 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 4406 #if defined(PETSC_HAVE_ELEMENTAL) 4407 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 4408 #endif 4409 #if defined(PETSC_HAVE_SCALAPACK) 4410 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*); 4411 #endif 4412 #if defined(PETSC_HAVE_HYPRE) 4413 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 4414 #endif 4415 4416 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 4417 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 4418 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat); 4419 4420 /*@C 4421 MatSeqAIJGetArray - gives read/write access to the array where the data for a MATSEQAIJ matrix is stored 4422 4423 Not Collective 4424 4425 Input Parameter: 4426 . mat - a MATSEQAIJ matrix 4427 4428 Output Parameter: 4429 . array - pointer to the data 4430 4431 Level: intermediate 4432 4433 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4434 @*/ 4435 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4436 { 4437 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4438 4439 PetscFunctionBegin; 4440 if (aij->ops->getarray) { 4441 PetscCall((*aij->ops->getarray)(A,array)); 4442 } else { 4443 *array = aij->a; 4444 } 4445 PetscFunctionReturn(0); 4446 } 4447 4448 /*@C 4449 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4450 4451 Not Collective 4452 4453 Input Parameters: 4454 + mat - a MATSEQAIJ matrix 4455 - array - pointer to the data 4456 4457 Level: intermediate 4458 4459 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayF90()` 4460 @*/ 4461 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4462 { 4463 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4464 4465 PetscFunctionBegin; 4466 if (aij->ops->restorearray) { 4467 PetscCall((*aij->ops->restorearray)(A,array)); 4468 } else { 4469 *array = NULL; 4470 } 4471 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4472 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4473 PetscFunctionReturn(0); 4474 } 4475 4476 /*@C 4477 MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a MATSEQAIJ matrix is stored 4478 4479 Not Collective 4480 4481 Input Parameter: 4482 . mat - a MATSEQAIJ matrix 4483 4484 Output Parameter: 4485 . array - pointer to the data 4486 4487 Level: intermediate 4488 4489 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4490 @*/ 4491 PetscErrorCode MatSeqAIJGetArrayRead(Mat A,const PetscScalar **array) 4492 { 4493 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4494 4495 PetscFunctionBegin; 4496 if (aij->ops->getarrayread) { 4497 PetscCall((*aij->ops->getarrayread)(A,array)); 4498 } else { 4499 *array = aij->a; 4500 } 4501 PetscFunctionReturn(0); 4502 } 4503 4504 /*@C 4505 MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4506 4507 Not Collective 4508 4509 Input Parameter: 4510 . mat - a MATSEQAIJ matrix 4511 4512 Output Parameter: 4513 . array - pointer to the data 4514 4515 Level: intermediate 4516 4517 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4518 @*/ 4519 PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A,const PetscScalar **array) 4520 { 4521 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4522 4523 PetscFunctionBegin; 4524 if (aij->ops->restorearrayread) { 4525 PetscCall((*aij->ops->restorearrayread)(A,array)); 4526 } else { 4527 *array = NULL; 4528 } 4529 PetscFunctionReturn(0); 4530 } 4531 4532 /*@C 4533 MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a MATSEQAIJ matrix is stored 4534 4535 Not Collective 4536 4537 Input Parameter: 4538 . mat - a MATSEQAIJ matrix 4539 4540 Output Parameter: 4541 . array - pointer to the data 4542 4543 Level: intermediate 4544 4545 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4546 @*/ 4547 PetscErrorCode MatSeqAIJGetArrayWrite(Mat A,PetscScalar **array) 4548 { 4549 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4550 4551 PetscFunctionBegin; 4552 if (aij->ops->getarraywrite) { 4553 PetscCall((*aij->ops->getarraywrite)(A,array)); 4554 } else { 4555 *array = aij->a; 4556 } 4557 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4558 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4559 PetscFunctionReturn(0); 4560 } 4561 4562 /*@C 4563 MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4564 4565 Not Collective 4566 4567 Input Parameter: 4568 . mat - a MATSEQAIJ matrix 4569 4570 Output Parameter: 4571 . array - pointer to the data 4572 4573 Level: intermediate 4574 4575 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4576 @*/ 4577 PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A,PetscScalar **array) 4578 { 4579 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4580 4581 PetscFunctionBegin; 4582 if (aij->ops->restorearraywrite) { 4583 PetscCall((*aij->ops->restorearraywrite)(A,array)); 4584 } else { 4585 *array = NULL; 4586 } 4587 PetscFunctionReturn(0); 4588 } 4589 4590 /*@C 4591 MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the SEQAIJ matrix 4592 4593 Not Collective 4594 4595 Input Parameter: 4596 . mat - a matrix of type MATSEQAIJ or its subclasses 4597 4598 Output Parameters: 4599 + i - row map array of the matrix 4600 . j - column index array of the matrix 4601 . a - data array of the matrix 4602 - memtype - memory type of the arrays 4603 4604 Notes: 4605 Any of the output parameters can be NULL, in which case the corresponding value is not returned. 4606 If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host. 4607 4608 One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix. 4609 If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix. 4610 4611 Level: Developer 4612 4613 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4614 @*/ 4615 PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat,const PetscInt **i,const PetscInt **j,PetscScalar **a,PetscMemType *mtype) 4616 { 4617 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 4618 4619 PetscFunctionBegin; 4620 PetscCheck(mat->preallocated,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"matrix is not preallocated"); 4621 if (aij->ops->getcsrandmemtype) { 4622 PetscCall((*aij->ops->getcsrandmemtype)(mat,i,j,a,mtype)); 4623 } else { 4624 if (i) *i = aij->i; 4625 if (j) *j = aij->j; 4626 if (a) *a = aij->a; 4627 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 4628 } 4629 PetscFunctionReturn(0); 4630 } 4631 4632 /*@C 4633 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4634 4635 Not Collective 4636 4637 Input Parameter: 4638 . mat - a MATSEQAIJ matrix 4639 4640 Output Parameter: 4641 . nz - the maximum number of nonzeros in any row 4642 4643 Level: intermediate 4644 4645 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4646 @*/ 4647 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4648 { 4649 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4650 4651 PetscFunctionBegin; 4652 *nz = aij->rmax; 4653 PetscFunctionReturn(0); 4654 } 4655 4656 PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, const PetscInt coo_i[], const PetscInt coo_j[]) 4657 { 4658 MPI_Comm comm; 4659 PetscInt *i,*j; 4660 PetscInt M,N,row; 4661 PetscCount k,p,q,nneg,nnz,start,end; /* Index the coo array, so use PetscCount as their type */ 4662 PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */ 4663 PetscInt *Aj; 4664 PetscScalar *Aa; 4665 Mat_SeqAIJ *seqaij = (Mat_SeqAIJ*)(mat->data); 4666 MatType rtype; 4667 PetscCount *perm,*jmap; 4668 4669 PetscFunctionBegin; 4670 PetscCall(MatResetPreallocationCOO_SeqAIJ(mat)); 4671 PetscCall(PetscObjectGetComm((PetscObject)mat,&comm)); 4672 PetscCall(MatGetSize(mat,&M,&N)); 4673 PetscCall(PetscMalloc2(coo_n,&i,coo_n,&j)); 4674 PetscCall(PetscArraycpy(i,coo_i,coo_n)); /* Make a copy since we'll modify it */ 4675 PetscCall(PetscArraycpy(j,coo_j,coo_n)); 4676 PetscCall(PetscMalloc1(coo_n,&perm)); 4677 for (k=0; k<coo_n; k++) { /* Ignore entries with negative row or col indices */ 4678 if (j[k] < 0) i[k] = -1; 4679 perm[k] = k; 4680 } 4681 4682 /* Sort by row */ 4683 PetscCall(PetscSortIntWithIntCountArrayPair(coo_n,i,j,perm)); 4684 for (k=0; k<coo_n; k++) {if (i[k] >= 0) break;} /* Advance k to the first row with a non-negative index */ 4685 nneg = k; 4686 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 */ 4687 nnz = 0; /* Total number of unique nonzeros to be counted */ 4688 jmap++; /* Inc jmap by 1 for convinience */ 4689 4690 PetscCall(PetscCalloc1(M+1,&Ai)); /* CSR of A */ 4691 PetscCall(PetscMalloc1(coo_n-nneg,&Aj)); /* We have at most coo_n-nneg unique nonzeros */ 4692 4693 /* In each row, sort by column, then unique column indices to get row length */ 4694 Ai++; /* Inc by 1 for convinience */ 4695 q = 0; /* q-th unique nonzero, with q starting from 0 */ 4696 while (k<coo_n) { 4697 row = i[k]; 4698 start = k; /* [start,end) indices for this row */ 4699 while (k<coo_n && i[k] == row) k++; 4700 end = k; 4701 PetscCall(PetscSortIntWithCountArray(end-start,j+start,perm+start)); 4702 /* Find number of unique col entries in this row */ 4703 Aj[q] = j[start]; /* Log the first nonzero in this row */ 4704 jmap[q] = 1; /* Number of repeats of this nozero entry */ 4705 Ai[row] = 1; 4706 nnz++; 4707 4708 for (p=start+1; p<end; p++) { /* Scan remaining nonzero in this row */ 4709 if (j[p] != j[p-1]) { /* Meet a new nonzero */ 4710 q++; 4711 jmap[q] = 1; 4712 Aj[q] = j[p]; 4713 Ai[row]++; 4714 nnz++; 4715 } else { 4716 jmap[q]++; 4717 } 4718 } 4719 q++; /* Move to next row and thus next unique nonzero */ 4720 } 4721 PetscCall(PetscFree2(i,j)); 4722 4723 Ai--; /* Back to the beginning of Ai[] */ 4724 for (k=0; k<M; k++) Ai[k+1] += Ai[k]; 4725 jmap--; /* Back to the beginning of jmap[] */ 4726 jmap[0] = 0; 4727 for (k=0; k<nnz; k++) jmap[k+1] += jmap[k]; 4728 if (nnz < coo_n-nneg) { /* Realloc with actual number of unique nonzeros */ 4729 PetscCount *jmap_new; 4730 PetscInt *Aj_new; 4731 4732 PetscCall(PetscMalloc1(nnz+1,&jmap_new)); 4733 PetscCall(PetscArraycpy(jmap_new,jmap,nnz+1)); 4734 PetscCall(PetscFree(jmap)); 4735 jmap = jmap_new; 4736 4737 PetscCall(PetscMalloc1(nnz,&Aj_new)); 4738 PetscCall(PetscArraycpy(Aj_new,Aj,nnz)); 4739 PetscCall(PetscFree(Aj)); 4740 Aj = Aj_new; 4741 } 4742 4743 if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */ 4744 PetscCount *perm_new; 4745 4746 PetscCall(PetscMalloc1(coo_n-nneg,&perm_new)); 4747 PetscCall(PetscArraycpy(perm_new,perm+nneg,coo_n-nneg)); 4748 PetscCall(PetscFree(perm)); 4749 perm = perm_new; 4750 } 4751 4752 PetscCall(MatGetRootType_Private(mat,&rtype)); 4753 PetscCall(PetscCalloc1(nnz,&Aa)); /* Zero the matrix */ 4754 PetscCall(MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF,M,N,Ai,Aj,Aa,rtype,mat)); 4755 4756 seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */ 4757 seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */ 4758 /* Record COO fields */ 4759 seqaij->coo_n = coo_n; 4760 seqaij->Atot = coo_n-nneg; /* Annz is seqaij->nz, so no need to record that again */ 4761 seqaij->jmap = jmap; /* of length nnz+1 */ 4762 seqaij->perm = perm; 4763 PetscFunctionReturn(0); 4764 } 4765 4766 static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A,const PetscScalar v[],InsertMode imode) 4767 { 4768 Mat_SeqAIJ *aseq = (Mat_SeqAIJ*)A->data; 4769 PetscCount i,j,Annz = aseq->nz; 4770 PetscCount *perm = aseq->perm,*jmap = aseq->jmap; 4771 PetscScalar *Aa; 4772 4773 PetscFunctionBegin; 4774 PetscCall(MatSeqAIJGetArray(A,&Aa)); 4775 for (i=0; i<Annz; i++) { 4776 PetscScalar sum = 0.0; 4777 for (j=jmap[i]; j<jmap[i+1]; j++) sum += v[perm[j]]; 4778 Aa[i] = (imode == INSERT_VALUES? 0.0 : Aa[i]) + sum; 4779 } 4780 PetscCall(MatSeqAIJRestoreArray(A,&Aa)); 4781 PetscFunctionReturn(0); 4782 } 4783 4784 #if defined(PETSC_HAVE_CUDA) 4785 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat,MatType,MatReuse,Mat*); 4786 #endif 4787 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4788 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat,MatType,MatReuse,Mat*); 4789 #endif 4790 4791 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4792 { 4793 Mat_SeqAIJ *b; 4794 PetscMPIInt size; 4795 4796 PetscFunctionBegin; 4797 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 4798 PetscCheck(size <= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4799 4800 PetscCall(PetscNewLog(B,&b)); 4801 4802 B->data = (void*)b; 4803 4804 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 4805 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4806 4807 b->row = NULL; 4808 b->col = NULL; 4809 b->icol = NULL; 4810 b->reallocs = 0; 4811 b->ignorezeroentries = PETSC_FALSE; 4812 b->roworiented = PETSC_TRUE; 4813 b->nonew = 0; 4814 b->diag = NULL; 4815 b->solve_work = NULL; 4816 B->spptr = NULL; 4817 b->saved_values = NULL; 4818 b->idiag = NULL; 4819 b->mdiag = NULL; 4820 b->ssor_work = NULL; 4821 b->omega = 1.0; 4822 b->fshift = 0.0; 4823 b->idiagvalid = PETSC_FALSE; 4824 b->ibdiagvalid = PETSC_FALSE; 4825 b->keepnonzeropattern = PETSC_FALSE; 4826 4827 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4828 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4829 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ)); 4830 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ)); 4831 #endif 4832 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ)); 4833 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ)); 4834 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ)); 4835 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ)); 4836 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ)); 4837 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM)); 4838 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL)); 4839 #if defined(PETSC_HAVE_MKL_SPARSE) 4840 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL)); 4841 #endif 4842 #if defined(PETSC_HAVE_CUDA) 4843 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 4844 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4845 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",MatProductSetFromOptions_SeqAIJ)); 4846 #endif 4847 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4848 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijkokkos_C",MatConvert_SeqAIJ_SeqAIJKokkos)); 4849 #endif 4850 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL)); 4851 #if defined(PETSC_HAVE_ELEMENTAL) 4852 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental)); 4853 #endif 4854 #if defined(PETSC_HAVE_SCALAPACK) 4855 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_scalapack_C",MatConvert_AIJ_ScaLAPACK)); 4856 #endif 4857 #if defined(PETSC_HAVE_HYPRE) 4858 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE)); 4859 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",MatProductSetFromOptions_Transpose_AIJ_AIJ)); 4860 #endif 4861 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense)); 4862 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL)); 4863 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS)); 4864 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ)); 4865 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ)); 4866 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ)); 4867 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ)); 4868 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ)); 4869 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ)); 4870 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_is_seqaij_C",MatProductSetFromOptions_IS_XAIJ)); 4871 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqaij_C",MatProductSetFromOptions_SeqDense_SeqAIJ)); 4872 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4873 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJKron_C",MatSeqAIJKron_SeqAIJ)); 4874 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_SeqAIJ)); 4875 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetValuesCOO_C",MatSetValuesCOO_SeqAIJ)); 4876 PetscCall(MatCreate_SeqAIJ_Inode(B)); 4877 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4878 PetscCall(MatSeqAIJSetTypeFromOptions(B)); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4879 PetscFunctionReturn(0); 4880 } 4881 4882 /* 4883 Given a matrix generated with MatGetFactor() duplicates all the information in A into C 4884 */ 4885 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4886 { 4887 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data,*a = (Mat_SeqAIJ*)A->data; 4888 PetscInt m = A->rmap->n,i; 4889 4890 PetscFunctionBegin; 4891 PetscCheck(A->assembled || cpvalues == MAT_DO_NOT_COPY_VALUES,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot duplicate unassembled matrix"); 4892 4893 C->factortype = A->factortype; 4894 c->row = NULL; 4895 c->col = NULL; 4896 c->icol = NULL; 4897 c->reallocs = 0; 4898 4899 C->assembled = A->assembled; 4900 C->preallocated = A->preallocated; 4901 4902 if (A->preallocated) { 4903 PetscCall(PetscLayoutReference(A->rmap,&C->rmap)); 4904 PetscCall(PetscLayoutReference(A->cmap,&C->cmap)); 4905 4906 PetscCall(PetscMalloc1(m,&c->imax)); 4907 PetscCall(PetscMemcpy(c->imax,a->imax,m*sizeof(PetscInt))); 4908 PetscCall(PetscMalloc1(m,&c->ilen)); 4909 PetscCall(PetscMemcpy(c->ilen,a->ilen,m*sizeof(PetscInt))); 4910 PetscCall(PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt))); 4911 4912 /* allocate the matrix space */ 4913 if (mallocmatspace) { 4914 PetscCall(PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i)); 4915 PetscCall(PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt))); 4916 4917 c->singlemalloc = PETSC_TRUE; 4918 4919 PetscCall(PetscArraycpy(c->i,a->i,m+1)); 4920 if (m > 0) { 4921 PetscCall(PetscArraycpy(c->j,a->j,a->i[m])); 4922 if (cpvalues == MAT_COPY_VALUES) { 4923 const PetscScalar *aa; 4924 4925 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4926 PetscCall(PetscArraycpy(c->a,aa,a->i[m])); 4927 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4928 } else { 4929 PetscCall(PetscArrayzero(c->a,a->i[m])); 4930 } 4931 } 4932 } 4933 4934 c->ignorezeroentries = a->ignorezeroentries; 4935 c->roworiented = a->roworiented; 4936 c->nonew = a->nonew; 4937 if (a->diag) { 4938 PetscCall(PetscMalloc1(m+1,&c->diag)); 4939 PetscCall(PetscMemcpy(c->diag,a->diag,m*sizeof(PetscInt))); 4940 PetscCall(PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt))); 4941 } else c->diag = NULL; 4942 4943 c->solve_work = NULL; 4944 c->saved_values = NULL; 4945 c->idiag = NULL; 4946 c->ssor_work = NULL; 4947 c->keepnonzeropattern = a->keepnonzeropattern; 4948 c->free_a = PETSC_TRUE; 4949 c->free_ij = PETSC_TRUE; 4950 4951 c->rmax = a->rmax; 4952 c->nz = a->nz; 4953 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4954 4955 c->compressedrow.use = a->compressedrow.use; 4956 c->compressedrow.nrows = a->compressedrow.nrows; 4957 if (a->compressedrow.use) { 4958 i = a->compressedrow.nrows; 4959 PetscCall(PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex)); 4960 PetscCall(PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1)); 4961 PetscCall(PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i)); 4962 } else { 4963 c->compressedrow.use = PETSC_FALSE; 4964 c->compressedrow.i = NULL; 4965 c->compressedrow.rindex = NULL; 4966 } 4967 c->nonzerorowcnt = a->nonzerorowcnt; 4968 C->nonzerostate = A->nonzerostate; 4969 4970 PetscCall(MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C)); 4971 } 4972 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist)); 4973 PetscFunctionReturn(0); 4974 } 4975 4976 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4977 { 4978 PetscFunctionBegin; 4979 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 4980 PetscCall(MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n)); 4981 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4982 PetscCall(MatSetBlockSizesFromMats(*B,A,A)); 4983 } 4984 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 4985 PetscCall(MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE)); 4986 PetscFunctionReturn(0); 4987 } 4988 4989 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4990 { 4991 PetscBool isbinary, ishdf5; 4992 4993 PetscFunctionBegin; 4994 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 4995 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 4996 /* force binary viewer to load .info file if it has not yet done so */ 4997 PetscCall(PetscViewerSetUp(viewer)); 4998 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 4999 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5)); 5000 if (isbinary) { 5001 PetscCall(MatLoad_SeqAIJ_Binary(newMat,viewer)); 5002 } else if (ishdf5) { 5003 #if defined(PETSC_HAVE_HDF5) 5004 PetscCall(MatLoad_AIJ_HDF5(newMat,viewer)); 5005 #else 5006 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 5007 #endif 5008 } else { 5009 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); 5010 } 5011 PetscFunctionReturn(0); 5012 } 5013 5014 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 5015 { 5016 Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data; 5017 PetscInt header[4],*rowlens,M,N,nz,sum,rows,cols,i; 5018 5019 PetscFunctionBegin; 5020 PetscCall(PetscViewerSetUp(viewer)); 5021 5022 /* read in matrix header */ 5023 PetscCall(PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT)); 5024 PetscCheck(header[0] == MAT_FILE_CLASSID,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 5025 M = header[1]; N = header[2]; nz = header[3]; 5026 PetscCheck(M >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%" PetscInt_FMT ") in file is negative",M); 5027 PetscCheck(N >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%" PetscInt_FMT ") in file is negative",N); 5028 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 5029 5030 /* set block sizes from the viewer's .info file */ 5031 PetscCall(MatLoad_Binary_BlockSizes(mat,viewer)); 5032 /* set local and global sizes if not set already */ 5033 if (mat->rmap->n < 0) mat->rmap->n = M; 5034 if (mat->cmap->n < 0) mat->cmap->n = N; 5035 if (mat->rmap->N < 0) mat->rmap->N = M; 5036 if (mat->cmap->N < 0) mat->cmap->N = N; 5037 PetscCall(PetscLayoutSetUp(mat->rmap)); 5038 PetscCall(PetscLayoutSetUp(mat->cmap)); 5039 5040 /* check if the matrix sizes are correct */ 5041 PetscCall(MatGetSize(mat,&rows,&cols)); 5042 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); 5043 5044 /* read in row lengths */ 5045 PetscCall(PetscMalloc1(M,&rowlens)); 5046 PetscCall(PetscViewerBinaryRead(viewer,rowlens,M,NULL,PETSC_INT)); 5047 /* check if sum(rowlens) is same as nz */ 5048 sum = 0; for (i=0; i<M; i++) sum += rowlens[i]; 5049 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); 5050 /* preallocate and check sizes */ 5051 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat,0,rowlens)); 5052 PetscCall(MatGetSize(mat,&rows,&cols)); 5053 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); 5054 /* store row lengths */ 5055 PetscCall(PetscArraycpy(a->ilen,rowlens,M)); 5056 PetscCall(PetscFree(rowlens)); 5057 5058 /* fill in "i" row pointers */ 5059 a->i[0] = 0; for (i=0; i<M; i++) a->i[i+1] = a->i[i] + a->ilen[i]; 5060 /* read in "j" column indices */ 5061 PetscCall(PetscViewerBinaryRead(viewer,a->j,nz,NULL,PETSC_INT)); 5062 /* read in "a" nonzero values */ 5063 PetscCall(PetscViewerBinaryRead(viewer,a->a,nz,NULL,PETSC_SCALAR)); 5064 5065 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 5066 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 5067 PetscFunctionReturn(0); 5068 } 5069 5070 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 5071 { 5072 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 5073 const PetscScalar *aa,*ba; 5074 #if defined(PETSC_USE_COMPLEX) 5075 PetscInt k; 5076 #endif 5077 5078 PetscFunctionBegin; 5079 /* If the matrix dimensions are not equal,or no of nonzeros */ 5080 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 5081 *flg = PETSC_FALSE; 5082 PetscFunctionReturn(0); 5083 } 5084 5085 /* if the a->i are the same */ 5086 PetscCall(PetscArraycmp(a->i,b->i,A->rmap->n+1,flg)); 5087 if (!*flg) PetscFunctionReturn(0); 5088 5089 /* if a->j are the same */ 5090 PetscCall(PetscArraycmp(a->j,b->j,a->nz,flg)); 5091 if (!*flg) PetscFunctionReturn(0); 5092 5093 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 5094 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 5095 /* if a->a are the same */ 5096 #if defined(PETSC_USE_COMPLEX) 5097 for (k=0; k<a->nz; k++) { 5098 if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) { 5099 *flg = PETSC_FALSE; 5100 PetscFunctionReturn(0); 5101 } 5102 } 5103 #else 5104 PetscCall(PetscArraycmp(aa,ba,a->nz,flg)); 5105 #endif 5106 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 5107 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 5108 PetscFunctionReturn(0); 5109 } 5110 5111 /*@ 5112 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 5113 provided by the user. 5114 5115 Collective 5116 5117 Input Parameters: 5118 + comm - must be an MPI communicator of size 1 5119 . m - number of rows 5120 . n - number of columns 5121 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 5122 . j - column indices 5123 - a - matrix values 5124 5125 Output Parameter: 5126 . mat - the matrix 5127 5128 Level: intermediate 5129 5130 Notes: 5131 The i, j, and a arrays are not copied by this routine, the user must free these arrays 5132 once the matrix is destroyed and not before 5133 5134 You cannot set new nonzero locations into this matrix, that will generate an error. 5135 5136 The i and j indices are 0 based 5137 5138 The format which is used for the sparse matrix input, is equivalent to a 5139 row-major ordering.. i.e for the following matrix, the input data expected is 5140 as shown 5141 5142 $ 1 0 0 5143 $ 2 0 3 5144 $ 4 5 6 5145 $ 5146 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 5147 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 5148 $ v = {1,2,3,4,5,6} [size = 6] 5149 5150 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateMPIAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()` 5151 5152 @*/ 5153 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 5154 { 5155 PetscInt ii; 5156 Mat_SeqAIJ *aij; 5157 PetscInt jj; 5158 5159 PetscFunctionBegin; 5160 PetscCheck(m <= 0 || i[0] == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 5161 PetscCall(MatCreate(comm,mat)); 5162 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5163 /* PetscCall(MatSetBlockSizes(*mat,,)); */ 5164 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5165 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,NULL)); 5166 aij = (Mat_SeqAIJ*)(*mat)->data; 5167 PetscCall(PetscMalloc1(m,&aij->imax)); 5168 PetscCall(PetscMalloc1(m,&aij->ilen)); 5169 5170 aij->i = i; 5171 aij->j = j; 5172 aij->a = a; 5173 aij->singlemalloc = PETSC_FALSE; 5174 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 5175 aij->free_a = PETSC_FALSE; 5176 aij->free_ij = PETSC_FALSE; 5177 5178 for (ii=0,aij->nonzerorowcnt=0,aij->rmax=0; ii<m; ii++) { 5179 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 5180 if (PetscDefined(USE_DEBUG)) { 5181 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]); 5182 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 5183 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); 5184 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); 5185 } 5186 } 5187 } 5188 if (PetscDefined(USE_DEBUG)) { 5189 for (ii=0; ii<aij->i[m]; ii++) { 5190 PetscCheck(j[ii] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 5191 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]); 5192 } 5193 } 5194 5195 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5196 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5197 PetscFunctionReturn(0); 5198 } 5199 5200 /*@ 5201 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 5202 provided by the user. 5203 5204 Collective 5205 5206 Input Parameters: 5207 + comm - must be an MPI communicator of size 1 5208 . m - number of rows 5209 . n - number of columns 5210 . i - row indices 5211 . j - column indices 5212 . a - matrix values 5213 . nz - number of nonzeros 5214 - idx - if the i and j indices start with 1 use PETSC_TRUE otherwise use PETSC_FALSE 5215 5216 Output Parameter: 5217 . mat - the matrix 5218 5219 Level: intermediate 5220 5221 Example: 5222 For the following matrix, the input data expected is as shown (using 0 based indexing) 5223 .vb 5224 1 0 0 5225 2 0 3 5226 4 5 6 5227 5228 i = {0,1,1,2,2,2} 5229 j = {0,0,2,0,1,2} 5230 v = {1,2,3,4,5,6} 5231 .ve 5232 5233 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateSeqAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()`, `MatSetValuesCOO()` 5234 5235 @*/ 5236 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 5237 { 5238 PetscInt ii, *nnz, one = 1,row,col; 5239 5240 PetscFunctionBegin; 5241 PetscCall(PetscCalloc1(m,&nnz)); 5242 for (ii = 0; ii < nz; ii++) { 5243 nnz[i[ii] - !!idx] += 1; 5244 } 5245 PetscCall(MatCreate(comm,mat)); 5246 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5247 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5248 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz)); 5249 for (ii = 0; ii < nz; ii++) { 5250 if (idx) { 5251 row = i[ii] - 1; 5252 col = j[ii] - 1; 5253 } else { 5254 row = i[ii]; 5255 col = j[ii]; 5256 } 5257 PetscCall(MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES)); 5258 } 5259 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5260 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5261 PetscCall(PetscFree(nnz)); 5262 PetscFunctionReturn(0); 5263 } 5264 5265 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 5266 { 5267 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 5268 5269 PetscFunctionBegin; 5270 a->idiagvalid = PETSC_FALSE; 5271 a->ibdiagvalid = PETSC_FALSE; 5272 5273 PetscCall(MatSeqAIJInvalidateDiagonal_Inode(A)); 5274 PetscFunctionReturn(0); 5275 } 5276 5277 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 5278 { 5279 PetscMPIInt size; 5280 5281 PetscFunctionBegin; 5282 PetscCallMPI(MPI_Comm_size(comm,&size)); 5283 if (size == 1) { 5284 if (scall == MAT_INITIAL_MATRIX) { 5285 PetscCall(MatDuplicate(inmat,MAT_COPY_VALUES,outmat)); 5286 } else { 5287 PetscCall(MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN)); 5288 } 5289 } else { 5290 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat)); 5291 } 5292 PetscFunctionReturn(0); 5293 } 5294 5295 /* 5296 Permute A into C's *local* index space using rowemb,colemb. 5297 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 5298 of [0,m), colemb is in [0,n). 5299 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 5300 */ 5301 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 5302 { 5303 /* If making this function public, change the error returned in this function away from _PLIB. */ 5304 Mat_SeqAIJ *Baij; 5305 PetscBool seqaij; 5306 PetscInt m,n,*nz,i,j,count; 5307 PetscScalar v; 5308 const PetscInt *rowindices,*colindices; 5309 5310 PetscFunctionBegin; 5311 if (!B) PetscFunctionReturn(0); 5312 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 5313 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij)); 5314 PetscCheck(seqaij,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 5315 if (rowemb) { 5316 PetscCall(ISGetLocalSize(rowemb,&m)); 5317 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); 5318 } else { 5319 PetscCheck(C->rmap->n == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 5320 } 5321 if (colemb) { 5322 PetscCall(ISGetLocalSize(colemb,&n)); 5323 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); 5324 } else { 5325 PetscCheck(C->cmap->n == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 5326 } 5327 5328 Baij = (Mat_SeqAIJ*)(B->data); 5329 if (pattern == DIFFERENT_NONZERO_PATTERN) { 5330 PetscCall(PetscMalloc1(B->rmap->n,&nz)); 5331 for (i=0; i<B->rmap->n; i++) { 5332 nz[i] = Baij->i[i+1] - Baij->i[i]; 5333 } 5334 PetscCall(MatSeqAIJSetPreallocation(C,0,nz)); 5335 PetscCall(PetscFree(nz)); 5336 } 5337 if (pattern == SUBSET_NONZERO_PATTERN) { 5338 PetscCall(MatZeroEntries(C)); 5339 } 5340 count = 0; 5341 rowindices = NULL; 5342 colindices = NULL; 5343 if (rowemb) { 5344 PetscCall(ISGetIndices(rowemb,&rowindices)); 5345 } 5346 if (colemb) { 5347 PetscCall(ISGetIndices(colemb,&colindices)); 5348 } 5349 for (i=0; i<B->rmap->n; i++) { 5350 PetscInt row; 5351 row = i; 5352 if (rowindices) row = rowindices[i]; 5353 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 5354 PetscInt col; 5355 col = Baij->j[count]; 5356 if (colindices) col = colindices[col]; 5357 v = Baij->a[count]; 5358 PetscCall(MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES)); 5359 ++count; 5360 } 5361 } 5362 /* FIXME: set C's nonzerostate correctly. */ 5363 /* Assembly for C is necessary. */ 5364 C->preallocated = PETSC_TRUE; 5365 C->assembled = PETSC_TRUE; 5366 C->was_assembled = PETSC_FALSE; 5367 PetscFunctionReturn(0); 5368 } 5369 5370 PetscFunctionList MatSeqAIJList = NULL; 5371 5372 /*@C 5373 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 5374 5375 Collective on Mat 5376 5377 Input Parameters: 5378 + mat - the matrix object 5379 - matype - matrix type 5380 5381 Options Database Key: 5382 . -mat_seqai_type <method> - for example seqaijcrl 5383 5384 Level: intermediate 5385 5386 .seealso: `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`, `Mat` 5387 @*/ 5388 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 5389 { 5390 PetscBool sametype; 5391 PetscErrorCode (*r)(Mat,MatType,MatReuse,Mat*); 5392 5393 PetscFunctionBegin; 5394 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 5395 PetscCall(PetscObjectTypeCompare((PetscObject)mat,matype,&sametype)); 5396 if (sametype) PetscFunctionReturn(0); 5397 5398 PetscCall(PetscFunctionListFind(MatSeqAIJList,matype,&r)); 5399 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 5400 PetscCall((*r)(mat,matype,MAT_INPLACE_MATRIX,&mat)); 5401 PetscFunctionReturn(0); 5402 } 5403 5404 /*@C 5405 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 5406 5407 Not Collective 5408 5409 Input Parameters: 5410 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 5411 - function - routine to convert to subtype 5412 5413 Notes: 5414 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 5415 5416 Then, your matrix can be chosen with the procedural interface at runtime via the option 5417 $ -mat_seqaij_type my_mat 5418 5419 Level: advanced 5420 5421 .seealso: `MatSeqAIJRegisterAll()` 5422 5423 Level: advanced 5424 @*/ 5425 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 5426 { 5427 PetscFunctionBegin; 5428 PetscCall(MatInitializePackage()); 5429 PetscCall(PetscFunctionListAdd(&MatSeqAIJList,sname,function)); 5430 PetscFunctionReturn(0); 5431 } 5432 5433 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 5434 5435 /*@C 5436 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 5437 5438 Not Collective 5439 5440 Level: advanced 5441 5442 .seealso: `MatRegisterAll()`, `MatSeqAIJRegister()` 5443 @*/ 5444 PetscErrorCode MatSeqAIJRegisterAll(void) 5445 { 5446 PetscFunctionBegin; 5447 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 5448 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 5449 5450 PetscCall(MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL)); 5451 PetscCall(MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM)); 5452 PetscCall(MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL)); 5453 #if defined(PETSC_HAVE_MKL_SPARSE) 5454 PetscCall(MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL)); 5455 #endif 5456 #if defined(PETSC_HAVE_CUDA) 5457 PetscCall(MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 5458 #endif 5459 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 5460 PetscCall(MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos)); 5461 #endif 5462 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 5463 PetscCall(MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL)); 5464 #endif 5465 PetscFunctionReturn(0); 5466 } 5467 5468 /* 5469 Special version for direct calls from Fortran 5470 */ 5471 #include <petsc/private/fortranimpl.h> 5472 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5473 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 5474 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 5475 #define matsetvaluesseqaij_ matsetvaluesseqaij 5476 #endif 5477 5478 /* Change these macros so can be used in void function */ 5479 5480 /* Change these macros so can be used in void function */ 5481 /* Identical to PetscCallVoid, except it assigns to *_ierr */ 5482 #undef PetscCall 5483 #define PetscCall(...) do { \ 5484 PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \ 5485 if (PetscUnlikely(ierr_msv_mpiaij)) { \ 5486 *_ierr = PetscError(PETSC_COMM_SELF,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr_msv_mpiaij,PETSC_ERROR_REPEAT," "); \ 5487 return; \ 5488 } \ 5489 } while (0) 5490 5491 #undef SETERRQ 5492 #define SETERRQ(comm,ierr,...) do { \ 5493 *_ierr = PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr,PETSC_ERROR_INITIAL,__VA_ARGS__); \ 5494 return; \ 5495 } while (0) 5496 5497 PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 5498 { 5499 Mat A = *AA; 5500 PetscInt m = *mm, n = *nn; 5501 InsertMode is = *isis; 5502 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 5503 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 5504 PetscInt *imax,*ai,*ailen; 5505 PetscInt *aj,nonew = a->nonew,lastcol = -1; 5506 MatScalar *ap,value,*aa; 5507 PetscBool ignorezeroentries = a->ignorezeroentries; 5508 PetscBool roworiented = a->roworiented; 5509 5510 PetscFunctionBegin; 5511 MatCheckPreallocated(A,1); 5512 imax = a->imax; 5513 ai = a->i; 5514 ailen = a->ilen; 5515 aj = a->j; 5516 aa = a->a; 5517 5518 for (k=0; k<m; k++) { /* loop over added rows */ 5519 row = im[k]; 5520 if (row < 0) continue; 5521 PetscCheck(row < A->rmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 5522 rp = aj + ai[row]; ap = aa + ai[row]; 5523 rmax = imax[row]; nrow = ailen[row]; 5524 low = 0; 5525 high = nrow; 5526 for (l=0; l<n; l++) { /* loop over added columns */ 5527 if (in[l] < 0) continue; 5528 PetscCheck(in[l] < A->cmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 5529 col = in[l]; 5530 if (roworiented) value = v[l + k*n]; 5531 else value = v[k + l*m]; 5532 5533 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 5534 5535 if (col <= lastcol) low = 0; 5536 else high = nrow; 5537 lastcol = col; 5538 while (high-low > 5) { 5539 t = (low+high)/2; 5540 if (rp[t] > col) high = t; 5541 else low = t; 5542 } 5543 for (i=low; i<high; i++) { 5544 if (rp[i] > col) break; 5545 if (rp[i] == col) { 5546 if (is == ADD_VALUES) ap[i] += value; 5547 else ap[i] = value; 5548 goto noinsert; 5549 } 5550 } 5551 if (value == 0.0 && ignorezeroentries) goto noinsert; 5552 if (nonew == 1) goto noinsert; 5553 PetscCheck(nonew != -1,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 5554 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 5555 N = nrow++ - 1; a->nz++; high++; 5556 /* shift up all the later entries in this row */ 5557 for (ii=N; ii>=i; ii--) { 5558 rp[ii+1] = rp[ii]; 5559 ap[ii+1] = ap[ii]; 5560 } 5561 rp[i] = col; 5562 ap[i] = value; 5563 A->nonzerostate++; 5564 noinsert:; 5565 low = i + 1; 5566 } 5567 ailen[row] = nrow; 5568 } 5569 PetscFunctionReturnVoid(); 5570 } 5571 /* Undefining these here since they were redefined from their original definition above! No 5572 * other PETSc functions should be defined past this point, as it is impossible to recover the 5573 * original definitions */ 5574 #undef PetscCall 5575 #undef SETERRQ 5576