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