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