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