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