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