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