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