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