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