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