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