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