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