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