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