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