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