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