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