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