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