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