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