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