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