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