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