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