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