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