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