1 #define PETSCMAT_DLL 2 3 /* 4 Defines the basic matrix operations for the AIJ (compressed row) 5 matrix storage format. 6 */ 7 8 #include "src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 9 #include "src/inline/spops.h" 10 #include "src/inline/dot.h" 11 #include "petscbt.h" 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 15 /* only works if matrix has a full set of diagonal entries */ 16 PetscErrorCode PETSCMAT_DLLEXPORT MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 17 { 18 PetscErrorCode ierr; 19 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 20 PetscInt i,*diag, m = Y->rmap.n; 21 PetscScalar *v,*aa = aij->a; 22 23 PetscFunctionBegin; 24 if (!Y->assembled) SETERRQ(PETSC_ERR_SUP,"Requires matrix be already assembled"); 25 ierr = MatMarkDiagonal_SeqAIJ(Y);CHKERRQ(ierr); 26 diag = aij->diag; 27 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 28 if (is == INSERT_VALUES) { 29 for (i=0; i<m; i++) { 30 aa[diag[i]] = v[i]; 31 } 32 } else { 33 for (i=0; i<m; i++) { 34 aa[diag[i]] += v[i]; 35 } 36 } 37 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 38 PetscFunctionReturn(0); 39 } 40 41 #undef __FUNCT__ 42 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 43 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 44 { 45 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 46 PetscErrorCode ierr; 47 PetscInt i,ishift; 48 49 PetscFunctionBegin; 50 *m = A->rmap.n; 51 if (!ia) PetscFunctionReturn(0); 52 ishift = 0; 53 if (symmetric && !A->structurally_symmetric) { 54 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap.n,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 55 } else if (oshift == 1) { 56 PetscInt nz = a->i[A->rmap.n]; 57 /* malloc space and add 1 to i and j indices */ 58 ierr = PetscMalloc((A->rmap.n+1)*sizeof(PetscInt),ia);CHKERRQ(ierr); 59 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),ja);CHKERRQ(ierr); 60 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 61 for (i=0; i<A->rmap.n+1; i++) (*ia)[i] = a->i[i] + 1; 62 } else { 63 *ia = a->i; *ja = a->j; 64 } 65 PetscFunctionReturn(0); 66 } 67 68 #undef __FUNCT__ 69 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 70 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 71 { 72 PetscErrorCode ierr; 73 74 PetscFunctionBegin; 75 if (!ia) PetscFunctionReturn(0); 76 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 77 ierr = PetscFree(*ia);CHKERRQ(ierr); 78 ierr = PetscFree(*ja);CHKERRQ(ierr); 79 } 80 PetscFunctionReturn(0); 81 } 82 83 #undef __FUNCT__ 84 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 85 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 86 { 87 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 88 PetscErrorCode ierr; 89 PetscInt i,*collengths,*cia,*cja,n = A->cmap.n,m = A->rmap.n; 90 PetscInt nz = a->i[m],row,*jj,mr,col; 91 92 PetscFunctionBegin; 93 *nn = n; 94 if (!ia) PetscFunctionReturn(0); 95 if (symmetric) { 96 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap.n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 97 } else { 98 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 99 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 100 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 101 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 102 jj = a->j; 103 for (i=0; i<nz; i++) { 104 collengths[jj[i]]++; 105 } 106 cia[0] = oshift; 107 for (i=0; i<n; i++) { 108 cia[i+1] = cia[i] + collengths[i]; 109 } 110 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 111 jj = a->j; 112 for (row=0; row<m; row++) { 113 mr = a->i[row+1] - a->i[row]; 114 for (i=0; i<mr; i++) { 115 col = *jj++; 116 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 117 } 118 } 119 ierr = PetscFree(collengths);CHKERRQ(ierr); 120 *ia = cia; *ja = cja; 121 } 122 PetscFunctionReturn(0); 123 } 124 125 #undef __FUNCT__ 126 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 127 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 128 { 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 if (!ia) PetscFunctionReturn(0); 133 134 ierr = PetscFree(*ia);CHKERRQ(ierr); 135 ierr = PetscFree(*ja);CHKERRQ(ierr); 136 137 PetscFunctionReturn(0); 138 } 139 140 #undef __FUNCT__ 141 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 142 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 143 { 144 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 145 PetscInt *ai = a->i; 146 PetscErrorCode ierr; 147 148 PetscFunctionBegin; 149 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 153 #define CHUNKSIZE 15 154 155 #undef __FUNCT__ 156 #define __FUNCT__ "MatSetValues_SeqAIJ" 157 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 158 { 159 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 160 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 161 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 162 PetscErrorCode ierr; 163 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 164 PetscScalar *ap,value,*aa = a->a; 165 PetscTruth ignorezeroentries = ((a->ignorezeroentries && is == ADD_VALUES) ? PETSC_TRUE:PETSC_FALSE); 166 PetscTruth roworiented = a->roworiented; 167 168 PetscFunctionBegin; 169 for (k=0; k<m; k++) { /* loop over added rows */ 170 row = im[k]; 171 if (row < 0) continue; 172 #if defined(PETSC_USE_DEBUG) 173 if (row >= A->rmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap.n-1); 174 #endif 175 rp = aj + ai[row]; ap = aa + ai[row]; 176 rmax = imax[row]; nrow = ailen[row]; 177 low = 0; 178 high = nrow; 179 for (l=0; l<n; l++) { /* loop over added columns */ 180 if (in[l] < 0) continue; 181 #if defined(PETSC_USE_DEBUG) 182 if (in[l] >= A->cmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap.n-1); 183 #endif 184 col = in[l]; 185 if (roworiented) { 186 value = v[l + k*n]; 187 } else { 188 value = v[k + l*m]; 189 } 190 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 191 192 if (col <= lastcol) low = 0; else high = nrow; 193 lastcol = col; 194 while (high-low > 5) { 195 t = (low+high)/2; 196 if (rp[t] > col) high = t; 197 else low = t; 198 } 199 for (i=low; i<high; i++) { 200 if (rp[i] > col) break; 201 if (rp[i] == col) { 202 if (is == ADD_VALUES) ap[i] += value; 203 else ap[i] = value; 204 goto noinsert; 205 } 206 } 207 if (value == 0.0 && ignorezeroentries) goto noinsert; 208 if (nonew == 1) goto noinsert; 209 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 210 MatSeqXAIJReallocateAIJ(a,1,nrow,row,col,rmax,aa,ai,aj,A->rmap.n,rp,ap,imax,nonew); 211 N = nrow++ - 1; a->nz++; high++; 212 /* shift up all the later entries in this row */ 213 for (ii=N; ii>=i; ii--) { 214 rp[ii+1] = rp[ii]; 215 ap[ii+1] = ap[ii]; 216 } 217 rp[i] = col; 218 ap[i] = value; 219 noinsert:; 220 low = i + 1; 221 } 222 ailen[row] = nrow; 223 } 224 A->same_nonzero = PETSC_FALSE; 225 PetscFunctionReturn(0); 226 } 227 228 229 #undef __FUNCT__ 230 #define __FUNCT__ "MatGetValues_SeqAIJ" 231 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 232 { 233 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 234 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 235 PetscInt *ai = a->i,*ailen = a->ilen; 236 PetscScalar *ap,*aa = a->a,zero = 0.0; 237 238 PetscFunctionBegin; 239 for (k=0; k<m; k++) { /* loop over rows */ 240 row = im[k]; 241 if (row < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); 242 if (row >= A->rmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap.n-1); 243 rp = aj + ai[row]; ap = aa + ai[row]; 244 nrow = ailen[row]; 245 for (l=0; l<n; l++) { /* loop over columns */ 246 if (in[l] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); 247 if (in[l] >= A->cmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap.n-1); 248 col = in[l] ; 249 high = nrow; low = 0; /* assume unsorted */ 250 while (high-low > 5) { 251 t = (low+high)/2; 252 if (rp[t] > col) high = t; 253 else low = t; 254 } 255 for (i=low; i<high; i++) { 256 if (rp[i] > col) break; 257 if (rp[i] == col) { 258 *v++ = ap[i]; 259 goto finished; 260 } 261 } 262 *v++ = zero; 263 finished:; 264 } 265 } 266 PetscFunctionReturn(0); 267 } 268 269 270 #undef __FUNCT__ 271 #define __FUNCT__ "MatView_SeqAIJ_Binary" 272 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 273 { 274 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 275 PetscErrorCode ierr; 276 PetscInt i,*col_lens; 277 int fd; 278 279 PetscFunctionBegin; 280 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 281 ierr = PetscMalloc((4+A->rmap.n)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 282 col_lens[0] = MAT_FILE_COOKIE; 283 col_lens[1] = A->rmap.n; 284 col_lens[2] = A->cmap.n; 285 col_lens[3] = a->nz; 286 287 /* store lengths of each row and write (including header) to file */ 288 for (i=0; i<A->rmap.n; i++) { 289 col_lens[4+i] = a->i[i+1] - a->i[i]; 290 } 291 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap.n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 292 ierr = PetscFree(col_lens);CHKERRQ(ierr); 293 294 /* store column indices (zero start index) */ 295 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 296 297 /* store nonzero values */ 298 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 299 PetscFunctionReturn(0); 300 } 301 302 EXTERN PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 303 304 #undef __FUNCT__ 305 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 306 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 307 { 308 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 309 PetscErrorCode ierr; 310 PetscInt i,j,m = A->rmap.n,shift=0; 311 const char *name; 312 PetscViewerFormat format; 313 314 PetscFunctionBegin; 315 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 316 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 317 if (format == PETSC_VIEWER_ASCII_MATLAB) { 318 PetscInt nofinalvalue = 0; 319 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap.n-!shift)) { 320 nofinalvalue = 1; 321 } 322 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 323 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap.n);CHKERRQ(ierr); 324 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 325 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 326 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 327 328 for (i=0; i<m; i++) { 329 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 330 #if defined(PETSC_USE_COMPLEX) 331 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e + %18.16ei \n",i+1,a->j[j]+!shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 332 #else 333 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 334 #endif 335 } 336 } 337 if (nofinalvalue) { 338 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap.n,0.0);CHKERRQ(ierr); 339 } 340 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 341 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 342 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 343 PetscFunctionReturn(0); 344 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 345 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 346 for (i=0; i<m; i++) { 347 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 348 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 349 #if defined(PETSC_USE_COMPLEX) 350 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 351 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 352 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 353 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 354 } else if (PetscRealPart(a->a[j]) != 0.0) { 355 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 356 } 357 #else 358 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 359 #endif 360 } 361 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 362 } 363 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 364 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 365 PetscInt nzd=0,fshift=1,*sptr; 366 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 367 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 368 for (i=0; i<m; i++) { 369 sptr[i] = nzd+1; 370 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 371 if (a->j[j] >= i) { 372 #if defined(PETSC_USE_COMPLEX) 373 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 374 #else 375 if (a->a[j] != 0.0) nzd++; 376 #endif 377 } 378 } 379 } 380 sptr[m] = nzd+1; 381 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 382 for (i=0; i<m+1; i+=6) { 383 if (i+4<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5]);CHKERRQ(ierr);} 384 else if (i+3<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]);CHKERRQ(ierr);} 385 else if (i+2<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);CHKERRQ(ierr);} 386 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 387 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 388 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 389 } 390 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 391 ierr = PetscFree(sptr);CHKERRQ(ierr); 392 for (i=0; i<m; i++) { 393 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 394 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 395 } 396 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 397 } 398 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 399 for (i=0; i<m; i++) { 400 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 401 if (a->j[j] >= i) { 402 #if defined(PETSC_USE_COMPLEX) 403 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 404 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 405 } 406 #else 407 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 408 #endif 409 } 410 } 411 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 412 } 413 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 414 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 415 PetscInt cnt = 0,jcnt; 416 PetscScalar value; 417 418 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 419 for (i=0; i<m; i++) { 420 jcnt = 0; 421 for (j=0; j<A->cmap.n; j++) { 422 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 423 value = a->a[cnt++]; 424 jcnt++; 425 } else { 426 value = 0.0; 427 } 428 #if defined(PETSC_USE_COMPLEX) 429 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 430 #else 431 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 432 #endif 433 } 434 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 435 } 436 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 437 } else { 438 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 439 for (i=0; i<m; i++) { 440 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 441 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 442 #if defined(PETSC_USE_COMPLEX) 443 if (PetscImaginaryPart(a->a[j]) > 0.0) { 444 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 445 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 446 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 447 } else { 448 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 449 } 450 #else 451 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 452 #endif 453 } 454 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 455 } 456 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 457 } 458 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 459 PetscFunctionReturn(0); 460 } 461 462 #undef __FUNCT__ 463 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 464 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 465 { 466 Mat A = (Mat) Aa; 467 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 468 PetscErrorCode ierr; 469 PetscInt i,j,m = A->rmap.n,color; 470 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 471 PetscViewer viewer; 472 PetscViewerFormat format; 473 474 PetscFunctionBegin; 475 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 476 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 477 478 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 479 /* loop over matrix elements drawing boxes */ 480 481 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 482 /* Blue for negative, Cyan for zero and Red for positive */ 483 color = PETSC_DRAW_BLUE; 484 for (i=0; i<m; i++) { 485 y_l = m - i - 1.0; y_r = y_l + 1.0; 486 for (j=a->i[i]; j<a->i[i+1]; j++) { 487 x_l = a->j[j] ; x_r = x_l + 1.0; 488 #if defined(PETSC_USE_COMPLEX) 489 if (PetscRealPart(a->a[j]) >= 0.) continue; 490 #else 491 if (a->a[j] >= 0.) continue; 492 #endif 493 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 494 } 495 } 496 color = PETSC_DRAW_CYAN; 497 for (i=0; i<m; i++) { 498 y_l = m - i - 1.0; y_r = y_l + 1.0; 499 for (j=a->i[i]; j<a->i[i+1]; j++) { 500 x_l = a->j[j]; x_r = x_l + 1.0; 501 if (a->a[j] != 0.) continue; 502 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 503 } 504 } 505 color = PETSC_DRAW_RED; 506 for (i=0; i<m; i++) { 507 y_l = m - i - 1.0; y_r = y_l + 1.0; 508 for (j=a->i[i]; j<a->i[i+1]; j++) { 509 x_l = a->j[j]; x_r = x_l + 1.0; 510 #if defined(PETSC_USE_COMPLEX) 511 if (PetscRealPart(a->a[j]) <= 0.) continue; 512 #else 513 if (a->a[j] <= 0.) continue; 514 #endif 515 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 516 } 517 } 518 } else { 519 /* use contour shading to indicate magnitude of values */ 520 /* first determine max of all nonzero values */ 521 PetscInt nz = a->nz,count; 522 PetscDraw popup; 523 PetscReal scale; 524 525 for (i=0; i<nz; i++) { 526 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 527 } 528 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 529 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 530 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 531 count = 0; 532 for (i=0; i<m; i++) { 533 y_l = m - i - 1.0; y_r = y_l + 1.0; 534 for (j=a->i[i]; j<a->i[i+1]; j++) { 535 x_l = a->j[j]; x_r = x_l + 1.0; 536 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 537 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 538 count++; 539 } 540 } 541 } 542 PetscFunctionReturn(0); 543 } 544 545 #undef __FUNCT__ 546 #define __FUNCT__ "MatView_SeqAIJ_Draw" 547 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 548 { 549 PetscErrorCode ierr; 550 PetscDraw draw; 551 PetscReal xr,yr,xl,yl,h,w; 552 PetscTruth isnull; 553 554 PetscFunctionBegin; 555 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 556 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 557 if (isnull) PetscFunctionReturn(0); 558 559 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 560 xr = A->cmap.n; yr = A->rmap.n; h = yr/10.0; w = xr/10.0; 561 xr += w; yr += h; xl = -w; yl = -h; 562 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 563 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 564 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 565 PetscFunctionReturn(0); 566 } 567 568 #undef __FUNCT__ 569 #define __FUNCT__ "MatView_SeqAIJ" 570 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 571 { 572 PetscErrorCode ierr; 573 PetscTruth issocket,iascii,isbinary,isdraw; 574 575 PetscFunctionBegin; 576 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 577 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 578 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 579 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 580 if (iascii) { 581 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 582 #if defined(PETSC_USE_SOCKET_VIEWER) 583 } else if (issocket) { 584 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 585 ierr = PetscViewerSocketPutSparse_Private(viewer,A->rmap.n,A->cmap.n,a->nz,a->a,a->i,a->j);CHKERRQ(ierr); 586 #endif 587 } else if (isbinary) { 588 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 589 } else if (isdraw) { 590 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 591 } else { 592 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 593 } 594 ierr = MatView_Inode(A,viewer);CHKERRQ(ierr); 595 PetscFunctionReturn(0); 596 } 597 598 #undef __FUNCT__ 599 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 600 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 601 { 602 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 603 PetscErrorCode ierr; 604 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 605 PetscInt m = A->rmap.n,*ip,N,*ailen = a->ilen,rmax = 0; 606 PetscScalar *aa = a->a,*ap; 607 PetscReal ratio=0.6; 608 609 PetscFunctionBegin; 610 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 611 612 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 613 for (i=1; i<m; i++) { 614 /* move each row back by the amount of empty slots (fshift) before it*/ 615 fshift += imax[i-1] - ailen[i-1]; 616 rmax = PetscMax(rmax,ailen[i]); 617 if (fshift) { 618 ip = aj + ai[i] ; 619 ap = aa + ai[i] ; 620 N = ailen[i]; 621 for (j=0; j<N; j++) { 622 ip[j-fshift] = ip[j]; 623 ap[j-fshift] = ap[j]; 624 } 625 } 626 ai[i] = ai[i-1] + ailen[i-1]; 627 } 628 if (m) { 629 fshift += imax[m-1] - ailen[m-1]; 630 ai[m] = ai[m-1] + ailen[m-1]; 631 } 632 /* reset ilen and imax for each row */ 633 for (i=0; i<m; i++) { 634 ailen[i] = imax[i] = ai[i+1] - ai[i]; 635 } 636 a->nz = ai[m]; 637 638 /* diagonals may have moved, so kill the diagonal pointers */ 639 if (fshift && a->diag) { 640 ierr = PetscFree(a->diag);CHKERRQ(ierr); 641 ierr = PetscLogObjectMemory(A,-(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 642 a->diag = 0; 643 } 644 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap.n,fshift,a->nz);CHKERRQ(ierr); 645 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 646 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 647 a->reallocs = 0; 648 A->info.nz_unneeded = (double)fshift; 649 a->rmax = rmax; 650 651 /* check for zero rows. If found a large number of zero rows, use CompressedRow functions */ 652 ierr = Mat_CheckCompressedRow(A,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 653 A->same_nonzero = PETSC_TRUE; 654 655 ierr = MatAssemblyEnd_Inode(A,mode);CHKERRQ(ierr); 656 PetscFunctionReturn(0); 657 } 658 659 #undef __FUNCT__ 660 #define __FUNCT__ "MatRealPart_SeqAIJ" 661 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 662 { 663 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 664 PetscInt i,nz = a->nz; 665 PetscScalar *aa = a->a; 666 667 PetscFunctionBegin; 668 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 669 PetscFunctionReturn(0); 670 } 671 672 #undef __FUNCT__ 673 #define __FUNCT__ "MatImaginaryPart_SeqAIJ" 674 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 675 { 676 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 677 PetscInt i,nz = a->nz; 678 PetscScalar *aa = a->a; 679 680 PetscFunctionBegin; 681 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 682 PetscFunctionReturn(0); 683 } 684 685 #undef __FUNCT__ 686 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 687 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 688 { 689 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 690 PetscErrorCode ierr; 691 692 PetscFunctionBegin; 693 ierr = PetscMemzero(a->a,(a->i[A->rmap.n])*sizeof(PetscScalar));CHKERRQ(ierr); 694 PetscFunctionReturn(0); 695 } 696 697 #undef __FUNCT__ 698 #define __FUNCT__ "MatDestroy_SeqAIJ" 699 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 700 { 701 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 702 PetscErrorCode ierr; 703 704 PetscFunctionBegin; 705 #if defined(PETSC_USE_LOG) 706 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap.n,A->cmap.n,a->nz); 707 #endif 708 if (a->freedata){ 709 ierr = MatSeqXAIJFreeAIJ(a->singlemalloc,&a->a,&a->j,&a->i);CHKERRQ(ierr); 710 } 711 if (a->row) { 712 ierr = ISDestroy(a->row);CHKERRQ(ierr); 713 } 714 if (a->col) { 715 ierr = ISDestroy(a->col);CHKERRQ(ierr); 716 } 717 ierr = PetscFree(a->diag);CHKERRQ(ierr); 718 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 719 ierr = PetscFree(a->idiag);CHKERRQ(ierr); 720 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 721 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 722 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 723 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 724 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 725 if (a->compressedrow.use){ierr = PetscFree(a->compressedrow.i);} 726 727 ierr = MatDestroy_Inode(A);CHKERRQ(ierr); 728 729 ierr = PetscFree(a);CHKERRQ(ierr); 730 731 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 732 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 733 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 734 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqsbaij_C","",PETSC_NULL);CHKERRQ(ierr); 735 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 736 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatIsTranspose_C","",PETSC_NULL);CHKERRQ(ierr); 737 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 738 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C","",PETSC_NULL);CHKERRQ(ierr); 739 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatReorderForNonzeroDiagonal_C","",PETSC_NULL);CHKERRQ(ierr); 740 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 741 PetscFunctionReturn(0); 742 } 743 744 #undef __FUNCT__ 745 #define __FUNCT__ "MatCompress_SeqAIJ" 746 PetscErrorCode MatCompress_SeqAIJ(Mat A) 747 { 748 PetscFunctionBegin; 749 PetscFunctionReturn(0); 750 } 751 752 #undef __FUNCT__ 753 #define __FUNCT__ "MatSetOption_SeqAIJ" 754 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op) 755 { 756 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 757 PetscErrorCode ierr; 758 759 PetscFunctionBegin; 760 switch (op) { 761 case MAT_ROW_ORIENTED: 762 a->roworiented = PETSC_TRUE; 763 break; 764 case MAT_KEEP_ZEROED_ROWS: 765 a->keepzeroedrows = PETSC_TRUE; 766 break; 767 case MAT_COLUMN_ORIENTED: 768 a->roworiented = PETSC_FALSE; 769 break; 770 case MAT_COLUMNS_SORTED: 771 a->sorted = PETSC_TRUE; 772 break; 773 case MAT_COLUMNS_UNSORTED: 774 a->sorted = PETSC_FALSE; 775 break; 776 case MAT_NO_NEW_NONZERO_LOCATIONS: 777 a->nonew = 1; 778 break; 779 case MAT_NEW_NONZERO_LOCATION_ERR: 780 a->nonew = -1; 781 break; 782 case MAT_NEW_NONZERO_ALLOCATION_ERR: 783 a->nonew = -2; 784 break; 785 case MAT_YES_NEW_NONZERO_LOCATIONS: 786 a->nonew = 0; 787 break; 788 case MAT_IGNORE_ZERO_ENTRIES: 789 a->ignorezeroentries = PETSC_TRUE; 790 break; 791 case MAT_USE_COMPRESSEDROW: 792 a->compressedrow.use = PETSC_TRUE; 793 break; 794 case MAT_DO_NOT_USE_COMPRESSEDROW: 795 a->compressedrow.use = PETSC_FALSE; 796 break; 797 case MAT_ROWS_SORTED: 798 case MAT_ROWS_UNSORTED: 799 case MAT_YES_NEW_DIAGONALS: 800 case MAT_IGNORE_OFF_PROC_ENTRIES: 801 case MAT_USE_HASH_TABLE: 802 ierr = PetscInfo(A,"Option ignored\n");CHKERRQ(ierr); 803 break; 804 case MAT_NO_NEW_DIAGONALS: 805 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 806 default: 807 break; 808 } 809 ierr = MatSetOption_Inode(A,op);CHKERRQ(ierr); 810 PetscFunctionReturn(0); 811 } 812 813 #undef __FUNCT__ 814 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 815 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 816 { 817 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 818 PetscErrorCode ierr; 819 PetscInt i,j,n; 820 PetscScalar *x,zero = 0.0; 821 822 PetscFunctionBegin; 823 ierr = VecSet(v,zero);CHKERRQ(ierr); 824 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 825 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 826 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 827 for (i=0; i<A->rmap.n; i++) { 828 for (j=a->i[i]; j<a->i[i+1]; j++) { 829 if (a->j[j] == i) { 830 x[i] = a->a[j]; 831 break; 832 } 833 } 834 } 835 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 836 PetscFunctionReturn(0); 837 } 838 839 #undef __FUNCT__ 840 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 841 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 842 { 843 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 844 PetscScalar *x,*y; 845 PetscErrorCode ierr; 846 PetscInt m = A->rmap.n; 847 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 848 PetscScalar *v,alpha; 849 PetscInt n,i,*idx,*ii,*ridx=PETSC_NULL; 850 Mat_CompressedRow cprow = a->compressedrow; 851 PetscTruth usecprow = cprow.use; 852 #endif 853 854 PetscFunctionBegin; 855 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 856 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 857 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 858 859 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 860 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 861 #else 862 if (usecprow){ 863 m = cprow.nrows; 864 ii = cprow.i; 865 ridx = cprow.rindex; 866 } else { 867 ii = a->i; 868 } 869 for (i=0; i<m; i++) { 870 idx = a->j + ii[i] ; 871 v = a->a + ii[i] ; 872 n = ii[i+1] - ii[i]; 873 if (usecprow){ 874 alpha = x[ridx[i]]; 875 } else { 876 alpha = x[i]; 877 } 878 while (n-->0) {y[*idx++] += alpha * *v++;} 879 } 880 #endif 881 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 882 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 883 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 884 PetscFunctionReturn(0); 885 } 886 887 #undef __FUNCT__ 888 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 889 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 890 { 891 PetscScalar zero = 0.0; 892 PetscErrorCode ierr; 893 894 PetscFunctionBegin; 895 ierr = VecSet(yy,zero);CHKERRQ(ierr); 896 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 897 PetscFunctionReturn(0); 898 } 899 900 901 #undef __FUNCT__ 902 #define __FUNCT__ "MatMult_SeqAIJ" 903 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 904 { 905 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 906 PetscScalar *x,*y,*aa; 907 PetscErrorCode ierr; 908 PetscInt m=A->rmap.n,*aj,*ii; 909 PetscInt n,i,j,*ridx=PETSC_NULL; 910 PetscScalar sum; 911 PetscTruth usecprow=a->compressedrow.use; 912 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 913 PetscInt jrow; 914 #endif 915 916 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 917 #pragma disjoint(*x,*y,*aa) 918 #endif 919 920 PetscFunctionBegin; 921 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 922 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 923 aj = a->j; 924 aa = a->a; 925 ii = a->i; 926 if (usecprow){ /* use compressed row format */ 927 m = a->compressedrow.nrows; 928 ii = a->compressedrow.i; 929 ridx = a->compressedrow.rindex; 930 for (i=0; i<m; i++){ 931 n = ii[i+1] - ii[i]; 932 aj = a->j + ii[i]; 933 aa = a->a + ii[i]; 934 sum = 0.0; 935 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 936 y[*ridx++] = sum; 937 } 938 } else { /* do not use compressed row format */ 939 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 940 fortranmultaij_(&m,x,ii,aj,aa,y); 941 #else 942 for (i=0; i<m; i++) { 943 jrow = ii[i]; 944 n = ii[i+1] - jrow; 945 sum = 0.0; 946 for (j=0; j<n; j++) { 947 sum += aa[jrow]*x[aj[jrow]]; jrow++; 948 } 949 y[i] = sum; 950 } 951 #endif 952 } 953 ierr = PetscLogFlops(2*a->nz - m);CHKERRQ(ierr); 954 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 955 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 956 PetscFunctionReturn(0); 957 } 958 959 #undef __FUNCT__ 960 #define __FUNCT__ "MatMultAdd_SeqAIJ" 961 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 962 { 963 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 964 PetscScalar *x,*y,*z,*aa; 965 PetscErrorCode ierr; 966 PetscInt m = A->rmap.n,*aj,*ii; 967 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 968 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 969 PetscScalar sum; 970 PetscTruth usecprow=a->compressedrow.use; 971 #endif 972 973 PetscFunctionBegin; 974 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 975 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 976 if (zz != yy) { 977 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 978 } else { 979 z = y; 980 } 981 982 aj = a->j; 983 aa = a->a; 984 ii = a->i; 985 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 986 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 987 #else 988 if (usecprow){ /* use compressed row format */ 989 if (zz != yy){ 990 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 991 } 992 m = a->compressedrow.nrows; 993 ii = a->compressedrow.i; 994 ridx = a->compressedrow.rindex; 995 for (i=0; i<m; i++){ 996 n = ii[i+1] - ii[i]; 997 aj = a->j + ii[i]; 998 aa = a->a + ii[i]; 999 sum = y[*ridx]; 1000 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 1001 z[*ridx++] = sum; 1002 } 1003 } else { /* do not use compressed row format */ 1004 for (i=0; i<m; i++) { 1005 jrow = ii[i]; 1006 n = ii[i+1] - jrow; 1007 sum = y[i]; 1008 for (j=0; j<n; j++) { 1009 sum += aa[jrow]*x[aj[jrow]]; jrow++; 1010 } 1011 z[i] = sum; 1012 } 1013 } 1014 #endif 1015 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 1016 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1017 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1018 if (zz != yy) { 1019 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1020 } 1021 PetscFunctionReturn(0); 1022 } 1023 1024 /* 1025 Adds diagonal pointers to sparse matrix structure. 1026 */ 1027 #undef __FUNCT__ 1028 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1029 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1030 { 1031 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1032 PetscErrorCode ierr; 1033 PetscInt i,j,*diag,m = A->rmap.n; 1034 1035 PetscFunctionBegin; 1036 if (a->diag) PetscFunctionReturn(0); 1037 1038 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&diag);CHKERRQ(ierr); 1039 ierr = PetscLogObjectMemory(A,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 1040 for (i=0; i<A->rmap.n; i++) { 1041 diag[i] = a->i[i+1]; 1042 for (j=a->i[i]; j<a->i[i+1]; j++) { 1043 if (a->j[j] == i) { 1044 diag[i] = j; 1045 break; 1046 } 1047 } 1048 } 1049 a->diag = diag; 1050 PetscFunctionReturn(0); 1051 } 1052 1053 /* 1054 Checks for missing diagonals 1055 */ 1056 #undef __FUNCT__ 1057 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1058 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A) 1059 { 1060 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1061 PetscErrorCode ierr; 1062 PetscInt *diag,*jj = a->j,i; 1063 1064 PetscFunctionBegin; 1065 if (A->rmap.n > 0 && !jj) SETERRQ(PETSC_ERR_PLIB,"Matrix has no entries therefor is missing diagonal"); 1066 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1067 diag = a->diag; 1068 for (i=0; i<A->rmap.n; i++) { 1069 if (jj[diag[i]] != i) { 1070 SETERRQ1(PETSC_ERR_PLIB,"Matrix is missing diagonal number %D",i); 1071 } 1072 } 1073 PetscFunctionReturn(0); 1074 } 1075 1076 #undef __FUNCT__ 1077 #define __FUNCT__ "MatRelax_SeqAIJ" 1078 PetscErrorCode MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1079 { 1080 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1081 PetscScalar *x,d,*xs,sum,*t,scale,*idiag=0,*mdiag; 1082 const PetscScalar *v = a->a, *b, *bs,*xb, *ts; 1083 PetscErrorCode ierr; 1084 PetscInt n = A->cmap.n,m = A->rmap.n,i; 1085 const PetscInt *idx,*diag; 1086 1087 PetscFunctionBegin; 1088 its = its*lits; 1089 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 1090 1091 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1092 diag = a->diag; 1093 if (!a->idiag) { 1094 ierr = PetscMalloc(3*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1095 a->ssor = a->idiag + m; 1096 mdiag = a->ssor + m; 1097 1098 v = a->a; 1099 1100 /* this is wrong when fshift omega changes each iteration */ 1101 if (omega == 1.0 && !fshift) { 1102 for (i=0; i<m; i++) { 1103 mdiag[i] = v[diag[i]]; 1104 a->idiag[i] = 1.0/v[diag[i]]; 1105 } 1106 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1107 } else { 1108 for (i=0; i<m; i++) { 1109 mdiag[i] = v[diag[i]]; 1110 a->idiag[i] = omega/(fshift + v[diag[i]]); 1111 } 1112 ierr = PetscLogFlops(2*m);CHKERRQ(ierr); 1113 } 1114 } 1115 t = a->ssor; 1116 idiag = a->idiag; 1117 mdiag = a->idiag + 2*m; 1118 1119 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1120 if (xx != bb) { 1121 ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1122 } else { 1123 b = x; 1124 } 1125 1126 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1127 xs = x; 1128 if (flag == SOR_APPLY_UPPER) { 1129 /* apply (U + D/omega) to the vector */ 1130 bs = b; 1131 for (i=0; i<m; i++) { 1132 d = fshift + a->a[diag[i]]; 1133 n = a->i[i+1] - diag[i] - 1; 1134 idx = a->j + diag[i] + 1; 1135 v = a->a + diag[i] + 1; 1136 sum = b[i]*d/omega; 1137 SPARSEDENSEDOT(sum,bs,v,idx,n); 1138 x[i] = sum; 1139 } 1140 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1141 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1142 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1143 PetscFunctionReturn(0); 1144 } 1145 1146 1147 /* Let A = L + U + D; where L is lower trianglar, 1148 U is upper triangular, E is diagonal; This routine applies 1149 1150 (L + E)^{-1} A (U + E)^{-1} 1151 1152 to a vector efficiently using Eisenstat's trick. This is for 1153 the case of SSOR preconditioner, so E is D/omega where omega 1154 is the relaxation factor. 1155 */ 1156 1157 if (flag == SOR_APPLY_LOWER) { 1158 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1159 } else if (flag & SOR_EISENSTAT) { 1160 /* Let A = L + U + D; where L is lower trianglar, 1161 U is upper triangular, E is diagonal; This routine applies 1162 1163 (L + E)^{-1} A (U + E)^{-1} 1164 1165 to a vector efficiently using Eisenstat's trick. This is for 1166 the case of SSOR preconditioner, so E is D/omega where omega 1167 is the relaxation factor. 1168 */ 1169 scale = (2.0/omega) - 1.0; 1170 1171 /* x = (E + U)^{-1} b */ 1172 for (i=m-1; i>=0; i--) { 1173 n = a->i[i+1] - diag[i] - 1; 1174 idx = a->j + diag[i] + 1; 1175 v = a->a + diag[i] + 1; 1176 sum = b[i]; 1177 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1178 x[i] = sum*idiag[i]; 1179 } 1180 1181 /* t = b - (2*E - D)x */ 1182 v = a->a; 1183 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1184 1185 /* t = (E + L)^{-1}t */ 1186 ts = t; 1187 diag = a->diag; 1188 for (i=0; i<m; i++) { 1189 n = diag[i] - a->i[i]; 1190 idx = a->j + a->i[i]; 1191 v = a->a + a->i[i]; 1192 sum = t[i]; 1193 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1194 t[i] = sum*idiag[i]; 1195 /* x = x + t */ 1196 x[i] += t[i]; 1197 } 1198 1199 ierr = PetscLogFlops(6*m-1 + 2*a->nz);CHKERRQ(ierr); 1200 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1201 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1202 PetscFunctionReturn(0); 1203 } 1204 if (flag & SOR_ZERO_INITIAL_GUESS) { 1205 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1206 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1207 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)b); 1208 #else 1209 for (i=0; i<m; i++) { 1210 n = diag[i] - a->i[i]; 1211 idx = a->j + a->i[i]; 1212 v = a->a + a->i[i]; 1213 sum = b[i]; 1214 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1215 x[i] = sum*idiag[i]; 1216 } 1217 #endif 1218 xb = x; 1219 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1220 } else xb = b; 1221 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1222 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1223 for (i=0; i<m; i++) { 1224 x[i] *= mdiag[i]; 1225 } 1226 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1227 } 1228 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1229 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1230 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)xb); 1231 #else 1232 for (i=m-1; i>=0; i--) { 1233 n = a->i[i+1] - diag[i] - 1; 1234 idx = a->j + diag[i] + 1; 1235 v = a->a + diag[i] + 1; 1236 sum = xb[i]; 1237 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1238 x[i] = sum*idiag[i]; 1239 } 1240 #endif 1241 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1242 } 1243 its--; 1244 } 1245 while (its--) { 1246 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1247 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1248 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1249 #else 1250 for (i=0; i<m; i++) { 1251 d = fshift + a->a[diag[i]]; 1252 n = a->i[i+1] - a->i[i]; 1253 idx = a->j + a->i[i]; 1254 v = a->a + a->i[i]; 1255 sum = b[i]; 1256 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1257 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1258 } 1259 #endif 1260 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1261 } 1262 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1263 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1264 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1265 #else 1266 for (i=m-1; i>=0; i--) { 1267 d = fshift + a->a[diag[i]]; 1268 n = a->i[i+1] - a->i[i]; 1269 idx = a->j + a->i[i]; 1270 v = a->a + a->i[i]; 1271 sum = b[i]; 1272 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1273 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1274 } 1275 #endif 1276 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1277 } 1278 } 1279 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1280 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1281 PetscFunctionReturn(0); 1282 } 1283 1284 #undef __FUNCT__ 1285 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1286 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1287 { 1288 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1289 1290 PetscFunctionBegin; 1291 info->rows_global = (double)A->rmap.n; 1292 info->columns_global = (double)A->cmap.n; 1293 info->rows_local = (double)A->rmap.n; 1294 info->columns_local = (double)A->cmap.n; 1295 info->block_size = 1.0; 1296 info->nz_allocated = (double)a->maxnz; 1297 info->nz_used = (double)a->nz; 1298 info->nz_unneeded = (double)(a->maxnz - a->nz); 1299 info->assemblies = (double)A->num_ass; 1300 info->mallocs = (double)a->reallocs; 1301 info->memory = A->mem; 1302 if (A->factor) { 1303 info->fill_ratio_given = A->info.fill_ratio_given; 1304 info->fill_ratio_needed = A->info.fill_ratio_needed; 1305 info->factor_mallocs = A->info.factor_mallocs; 1306 } else { 1307 info->fill_ratio_given = 0; 1308 info->fill_ratio_needed = 0; 1309 info->factor_mallocs = 0; 1310 } 1311 PetscFunctionReturn(0); 1312 } 1313 1314 #undef __FUNCT__ 1315 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1316 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1317 { 1318 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1319 PetscInt i,m = A->rmap.n - 1; 1320 PetscErrorCode ierr; 1321 1322 PetscFunctionBegin; 1323 if (a->keepzeroedrows) { 1324 for (i=0; i<N; i++) { 1325 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1326 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1327 } 1328 if (diag != 0.0) { 1329 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1330 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1331 for (i=0; i<N; i++) { 1332 a->a[a->diag[rows[i]]] = diag; 1333 } 1334 } 1335 A->same_nonzero = PETSC_TRUE; 1336 } else { 1337 if (diag != 0.0) { 1338 for (i=0; i<N; i++) { 1339 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1340 if (a->ilen[rows[i]] > 0) { 1341 a->ilen[rows[i]] = 1; 1342 a->a[a->i[rows[i]]] = diag; 1343 a->j[a->i[rows[i]]] = rows[i]; 1344 } else { /* in case row was completely empty */ 1345 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1346 } 1347 } 1348 } else { 1349 for (i=0; i<N; i++) { 1350 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1351 a->ilen[rows[i]] = 0; 1352 } 1353 } 1354 A->same_nonzero = PETSC_FALSE; 1355 } 1356 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1357 PetscFunctionReturn(0); 1358 } 1359 1360 #undef __FUNCT__ 1361 #define __FUNCT__ "MatGetRow_SeqAIJ" 1362 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1363 { 1364 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1365 PetscInt *itmp; 1366 1367 PetscFunctionBegin; 1368 if (row < 0 || row >= A->rmap.n) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1369 1370 *nz = a->i[row+1] - a->i[row]; 1371 if (v) *v = a->a + a->i[row]; 1372 if (idx) { 1373 itmp = a->j + a->i[row]; 1374 if (*nz) { 1375 *idx = itmp; 1376 } 1377 else *idx = 0; 1378 } 1379 PetscFunctionReturn(0); 1380 } 1381 1382 /* remove this function? */ 1383 #undef __FUNCT__ 1384 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1385 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1386 { 1387 PetscFunctionBegin; 1388 PetscFunctionReturn(0); 1389 } 1390 1391 #undef __FUNCT__ 1392 #define __FUNCT__ "MatNorm_SeqAIJ" 1393 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1394 { 1395 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1396 PetscScalar *v = a->a; 1397 PetscReal sum = 0.0; 1398 PetscErrorCode ierr; 1399 PetscInt i,j; 1400 1401 PetscFunctionBegin; 1402 if (type == NORM_FROBENIUS) { 1403 for (i=0; i<a->nz; i++) { 1404 #if defined(PETSC_USE_COMPLEX) 1405 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1406 #else 1407 sum += (*v)*(*v); v++; 1408 #endif 1409 } 1410 *nrm = sqrt(sum); 1411 } else if (type == NORM_1) { 1412 PetscReal *tmp; 1413 PetscInt *jj = a->j; 1414 ierr = PetscMalloc((A->cmap.n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1415 ierr = PetscMemzero(tmp,A->cmap.n*sizeof(PetscReal));CHKERRQ(ierr); 1416 *nrm = 0.0; 1417 for (j=0; j<a->nz; j++) { 1418 tmp[*jj++] += PetscAbsScalar(*v); v++; 1419 } 1420 for (j=0; j<A->cmap.n; j++) { 1421 if (tmp[j] > *nrm) *nrm = tmp[j]; 1422 } 1423 ierr = PetscFree(tmp);CHKERRQ(ierr); 1424 } else if (type == NORM_INFINITY) { 1425 *nrm = 0.0; 1426 for (j=0; j<A->rmap.n; j++) { 1427 v = a->a + a->i[j]; 1428 sum = 0.0; 1429 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1430 sum += PetscAbsScalar(*v); v++; 1431 } 1432 if (sum > *nrm) *nrm = sum; 1433 } 1434 } else { 1435 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1436 } 1437 PetscFunctionReturn(0); 1438 } 1439 1440 #undef __FUNCT__ 1441 #define __FUNCT__ "MatTranspose_SeqAIJ" 1442 PetscErrorCode MatTranspose_SeqAIJ(Mat A,Mat *B) 1443 { 1444 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1445 Mat C; 1446 PetscErrorCode ierr; 1447 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap.n,len,*col; 1448 PetscScalar *array = a->a; 1449 1450 PetscFunctionBegin; 1451 if (!B && m != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1452 ierr = PetscMalloc((1+A->cmap.n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1453 ierr = PetscMemzero(col,(1+A->cmap.n)*sizeof(PetscInt));CHKERRQ(ierr); 1454 1455 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1456 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1457 ierr = MatSetSizes(C,A->cmap.n,m,A->cmap.n,m);CHKERRQ(ierr); 1458 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1459 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1460 ierr = PetscFree(col);CHKERRQ(ierr); 1461 for (i=0; i<m; i++) { 1462 len = ai[i+1]-ai[i]; 1463 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1464 array += len; 1465 aj += len; 1466 } 1467 1468 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1469 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1470 1471 if (B) { 1472 *B = C; 1473 } else { 1474 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1475 } 1476 PetscFunctionReturn(0); 1477 } 1478 1479 EXTERN_C_BEGIN 1480 #undef __FUNCT__ 1481 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1482 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1483 { 1484 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1485 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1486 PetscErrorCode ierr; 1487 PetscInt ma,na,mb,nb, i; 1488 1489 PetscFunctionBegin; 1490 bij = (Mat_SeqAIJ *) B->data; 1491 1492 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1493 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1494 if (ma!=nb || na!=mb){ 1495 *f = PETSC_FALSE; 1496 PetscFunctionReturn(0); 1497 } 1498 aii = aij->i; bii = bij->i; 1499 adx = aij->j; bdx = bij->j; 1500 va = aij->a; vb = bij->a; 1501 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1502 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1503 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1504 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1505 1506 *f = PETSC_TRUE; 1507 for (i=0; i<ma; i++) { 1508 while (aptr[i]<aii[i+1]) { 1509 PetscInt idc,idr; 1510 PetscScalar vc,vr; 1511 /* column/row index/value */ 1512 idc = adx[aptr[i]]; 1513 idr = bdx[bptr[idc]]; 1514 vc = va[aptr[i]]; 1515 vr = vb[bptr[idc]]; 1516 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1517 *f = PETSC_FALSE; 1518 goto done; 1519 } else { 1520 aptr[i]++; 1521 if (B || i!=idc) bptr[idc]++; 1522 } 1523 } 1524 } 1525 done: 1526 ierr = PetscFree(aptr);CHKERRQ(ierr); 1527 if (B) { 1528 ierr = PetscFree(bptr);CHKERRQ(ierr); 1529 } 1530 PetscFunctionReturn(0); 1531 } 1532 EXTERN_C_END 1533 1534 #undef __FUNCT__ 1535 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1536 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1537 { 1538 PetscErrorCode ierr; 1539 PetscFunctionBegin; 1540 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1541 PetscFunctionReturn(0); 1542 } 1543 1544 #undef __FUNCT__ 1545 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1546 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1547 { 1548 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1549 PetscScalar *l,*r,x,*v; 1550 PetscErrorCode ierr; 1551 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,M,nz = a->nz,*jj; 1552 1553 PetscFunctionBegin; 1554 if (ll) { 1555 /* The local size is used so that VecMPI can be passed to this routine 1556 by MatDiagonalScale_MPIAIJ */ 1557 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1558 if (m != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1559 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1560 v = a->a; 1561 for (i=0; i<m; i++) { 1562 x = l[i]; 1563 M = a->i[i+1] - a->i[i]; 1564 for (j=0; j<M; j++) { (*v++) *= x;} 1565 } 1566 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1567 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1568 } 1569 if (rr) { 1570 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1571 if (n != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1572 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1573 v = a->a; jj = a->j; 1574 for (i=0; i<nz; i++) { 1575 (*v++) *= r[*jj++]; 1576 } 1577 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1578 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1579 } 1580 PetscFunctionReturn(0); 1581 } 1582 1583 #undef __FUNCT__ 1584 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1585 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1586 { 1587 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1588 PetscErrorCode ierr; 1589 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap.n,*lens; 1590 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1591 PetscInt *irow,*icol,nrows,ncols; 1592 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1593 PetscScalar *a_new,*mat_a; 1594 Mat C; 1595 PetscTruth stride; 1596 1597 PetscFunctionBegin; 1598 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1599 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1600 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1601 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1602 1603 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1604 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1605 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1606 1607 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1608 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1609 if (stride && step == 1) { 1610 /* special case of contiguous rows */ 1611 ierr = PetscMalloc((2*nrows+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1612 starts = lens + nrows; 1613 /* loop over new rows determining lens and starting points */ 1614 for (i=0; i<nrows; i++) { 1615 kstart = ai[irow[i]]; 1616 kend = kstart + ailen[irow[i]]; 1617 for (k=kstart; k<kend; k++) { 1618 if (aj[k] >= first) { 1619 starts[i] = k; 1620 break; 1621 } 1622 } 1623 sum = 0; 1624 while (k < kend) { 1625 if (aj[k++] >= first+ncols) break; 1626 sum++; 1627 } 1628 lens[i] = sum; 1629 } 1630 /* create submatrix */ 1631 if (scall == MAT_REUSE_MATRIX) { 1632 PetscInt n_cols,n_rows; 1633 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1634 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1635 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1636 C = *B; 1637 } else { 1638 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1639 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1640 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1641 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1642 } 1643 c = (Mat_SeqAIJ*)C->data; 1644 1645 /* loop over rows inserting into submatrix */ 1646 a_new = c->a; 1647 j_new = c->j; 1648 i_new = c->i; 1649 1650 for (i=0; i<nrows; i++) { 1651 ii = starts[i]; 1652 lensi = lens[i]; 1653 for (k=0; k<lensi; k++) { 1654 *j_new++ = aj[ii+k] - first; 1655 } 1656 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1657 a_new += lensi; 1658 i_new[i+1] = i_new[i] + lensi; 1659 c->ilen[i] = lensi; 1660 } 1661 ierr = PetscFree(lens);CHKERRQ(ierr); 1662 } else { 1663 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1664 ierr = PetscMalloc((1+oldcols)*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1665 1666 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1667 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1668 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1669 /* determine lens of each row */ 1670 for (i=0; i<nrows; i++) { 1671 kstart = ai[irow[i]]; 1672 kend = kstart + a->ilen[irow[i]]; 1673 lens[i] = 0; 1674 for (k=kstart; k<kend; k++) { 1675 if (smap[aj[k]]) { 1676 lens[i]++; 1677 } 1678 } 1679 } 1680 /* Create and fill new matrix */ 1681 if (scall == MAT_REUSE_MATRIX) { 1682 PetscTruth equal; 1683 1684 c = (Mat_SeqAIJ *)((*B)->data); 1685 if ((*B)->rmap.n != nrows || (*B)->cmap.n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1686 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap.n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1687 if (!equal) { 1688 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1689 } 1690 ierr = PetscMemzero(c->ilen,(*B)->rmap.n*sizeof(PetscInt));CHKERRQ(ierr); 1691 C = *B; 1692 } else { 1693 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1694 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1695 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1696 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1697 } 1698 c = (Mat_SeqAIJ *)(C->data); 1699 for (i=0; i<nrows; i++) { 1700 row = irow[i]; 1701 kstart = ai[row]; 1702 kend = kstart + a->ilen[row]; 1703 mat_i = c->i[i]; 1704 mat_j = c->j + mat_i; 1705 mat_a = c->a + mat_i; 1706 mat_ilen = c->ilen + i; 1707 for (k=kstart; k<kend; k++) { 1708 if ((tcol=smap[a->j[k]])) { 1709 *mat_j++ = tcol - 1; 1710 *mat_a++ = a->a[k]; 1711 (*mat_ilen)++; 1712 1713 } 1714 } 1715 } 1716 /* Free work space */ 1717 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1718 ierr = PetscFree(smap);CHKERRQ(ierr); 1719 ierr = PetscFree(lens);CHKERRQ(ierr); 1720 } 1721 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1722 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1723 1724 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1725 *B = C; 1726 PetscFunctionReturn(0); 1727 } 1728 1729 /* 1730 */ 1731 #undef __FUNCT__ 1732 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1733 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1734 { 1735 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1736 PetscErrorCode ierr; 1737 Mat outA; 1738 PetscTruth row_identity,col_identity; 1739 1740 PetscFunctionBegin; 1741 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1742 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1743 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1744 if (!row_identity || !col_identity) { 1745 SETERRQ(PETSC_ERR_ARG_WRONG,"Row and column permutations must be identity for in-place ILU"); 1746 } 1747 1748 outA = inA; 1749 inA->factor = FACTOR_LU; 1750 a->row = row; 1751 a->col = col; 1752 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1753 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1754 1755 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1756 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1757 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1758 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1759 1760 if (!a->solve_work) { /* this matrix may have been factored before */ 1761 ierr = PetscMalloc((inA->rmap.n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1762 } 1763 1764 if (!a->diag) { 1765 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1766 } 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*/ MatFDColoringApply_SeqAIJ, 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->freedata = PETSC_TRUE; 2744 } else { 2745 b->freedata = PETSC_FALSE; 2746 } 2747 2748 b->nz = 0; 2749 b->maxnz = nz; 2750 B->info.nz_unneeded = (double)b->maxnz; 2751 PetscFunctionReturn(0); 2752 } 2753 EXTERN_C_END 2754 2755 #undef __FUNCT__ 2756 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 2757 /*@C 2758 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 2759 2760 Input Parameters: 2761 + B - the matrix 2762 . i - the indices into j for the start of each row (starts with zero) 2763 . j - the column indices for each row (starts with zero) these must be sorted for each row 2764 - v - optional values in the matrix 2765 2766 Contributed by: Lisandro Dalchin 2767 2768 Level: developer 2769 2770 .keywords: matrix, aij, compressed row, sparse, sequential 2771 2772 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 2773 @*/ 2774 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 2775 { 2776 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 2777 PetscErrorCode ierr; 2778 2779 PetscFunctionBegin; 2780 PetscValidHeaderSpecific(B,MAT_COOKIE,1); 2781 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 2782 if (f) { 2783 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 2784 } 2785 PetscFunctionReturn(0); 2786 } 2787 2788 EXTERN_C_BEGIN 2789 #undef __FUNCT__ 2790 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 2791 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt I[],const PetscInt J[],const PetscScalar v[]) 2792 { 2793 PetscInt i; 2794 PetscInt m,n; 2795 PetscInt nz; 2796 PetscInt *nnz, nz_max = 0; 2797 PetscScalar *values; 2798 PetscErrorCode ierr; 2799 2800 PetscFunctionBegin; 2801 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 2802 2803 if (I[0]) { 2804 SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "I[0] must be 0 it is %D", I[0]); 2805 } 2806 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 2807 for(i = 0; i < m; i++) { 2808 nz = I[i+1]- I[i]; 2809 nz_max = PetscMax(nz_max, nz); 2810 if (nz < 0) { 2811 SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 2812 } 2813 nnz[i] = nz; 2814 } 2815 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 2816 ierr = PetscFree(nnz);CHKERRQ(ierr); 2817 2818 if (v) { 2819 values = (PetscScalar*) v; 2820 } else { 2821 ierr = PetscMalloc((nz_max+1)*sizeof(PetscScalar), &values);CHKERRQ(ierr); 2822 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 2823 } 2824 2825 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 2826 2827 for(i = 0; i < m; i++) { 2828 nz = I[i+1] - I[i]; 2829 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+I[i], values + (v ? I[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 2830 } 2831 2832 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2833 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2834 ierr = MatSetOption(B,MAT_COLUMNS_UNSORTED);CHKERRQ(ierr); 2835 2836 if (!v) { 2837 ierr = PetscFree(values);CHKERRQ(ierr); 2838 } 2839 PetscFunctionReturn(0); 2840 } 2841 EXTERN_C_END 2842 2843 /*MC 2844 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 2845 based on compressed sparse row format. 2846 2847 Options Database Keys: 2848 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 2849 2850 Level: beginner 2851 2852 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 2853 M*/ 2854 2855 EXTERN_C_BEGIN 2856 #undef __FUNCT__ 2857 #define __FUNCT__ "MatCreate_SeqAIJ" 2858 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 2859 { 2860 Mat_SeqAIJ *b; 2861 PetscErrorCode ierr; 2862 PetscMPIInt size; 2863 2864 PetscFunctionBegin; 2865 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2866 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2867 2868 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2869 B->data = (void*)b; 2870 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2871 B->factor = 0; 2872 B->mapping = 0; 2873 b->row = 0; 2874 b->col = 0; 2875 b->icol = 0; 2876 b->reallocs = 0; 2877 b->sorted = PETSC_FALSE; 2878 b->ignorezeroentries = PETSC_FALSE; 2879 b->roworiented = PETSC_TRUE; 2880 b->nonew = 0; 2881 b->diag = 0; 2882 b->solve_work = 0; 2883 B->spptr = 0; 2884 b->saved_values = 0; 2885 b->idiag = 0; 2886 b->ssor = 0; 2887 b->keepzeroedrows = PETSC_FALSE; 2888 b->xtoy = 0; 2889 b->XtoY = 0; 2890 b->compressedrow.use = PETSC_FALSE; 2891 b->compressedrow.nrows = B->rmap.n; 2892 b->compressedrow.i = PETSC_NULL; 2893 b->compressedrow.rindex = PETSC_NULL; 2894 b->compressedrow.checked = PETSC_FALSE; 2895 B->same_nonzero = PETSC_FALSE; 2896 2897 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2898 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2899 "MatSeqAIJSetColumnIndices_SeqAIJ", 2900 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2901 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2902 "MatStoreValues_SeqAIJ", 2903 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2904 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2905 "MatRetrieveValues_SeqAIJ", 2906 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2907 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2908 "MatConvert_SeqAIJ_SeqSBAIJ", 2909 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2910 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2911 "MatConvert_SeqAIJ_SeqBAIJ", 2912 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2913 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 2914 "MatIsTranspose_SeqAIJ", 2915 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 2916 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2917 "MatSeqAIJSetPreallocation_SeqAIJ", 2918 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2919 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 2920 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 2921 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 2922 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2923 "MatReorderForNonzeroDiagonal_SeqAIJ", 2924 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2925 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2926 "MatMatMult_SeqAIJ_SeqDense", 2927 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2928 ierr = MatCreate_Inode(B);CHKERRQ(ierr); 2929 PetscFunctionReturn(0); 2930 } 2931 EXTERN_C_END 2932 2933 #undef __FUNCT__ 2934 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2935 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2936 { 2937 Mat C; 2938 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2939 PetscErrorCode ierr; 2940 PetscInt i,m = A->rmap.n; 2941 2942 PetscFunctionBegin; 2943 *B = 0; 2944 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 2945 ierr = MatSetSizes(C,A->rmap.n,A->cmap.n,A->rmap.n,A->cmap.n);CHKERRQ(ierr); 2946 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 2947 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2948 2949 ierr = PetscMapCopy(A->comm,&A->rmap,&C->rmap);CHKERRQ(ierr); 2950 ierr = PetscMapCopy(A->comm,&A->cmap,&C->cmap);CHKERRQ(ierr); 2951 2952 c = (Mat_SeqAIJ*)C->data; 2953 2954 C->factor = A->factor; 2955 2956 c->row = 0; 2957 c->col = 0; 2958 c->icol = 0; 2959 c->reallocs = 0; 2960 2961 C->assembled = PETSC_TRUE; 2962 2963 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 2964 for (i=0; i<m; i++) { 2965 c->imax[i] = a->imax[i]; 2966 c->ilen[i] = a->ilen[i]; 2967 } 2968 2969 /* allocate the matrix space */ 2970 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 2971 c->singlemalloc = PETSC_TRUE; 2972 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2973 if (m > 0) { 2974 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 2975 if (cpvalues == MAT_COPY_VALUES) { 2976 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2977 } else { 2978 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2979 } 2980 } 2981 2982 c->sorted = a->sorted; 2983 c->ignorezeroentries = a->ignorezeroentries; 2984 c->roworiented = a->roworiented; 2985 c->nonew = a->nonew; 2986 if (a->diag) { 2987 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2988 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2989 for (i=0; i<m; i++) { 2990 c->diag[i] = a->diag[i]; 2991 } 2992 } else c->diag = 0; 2993 c->solve_work = 0; 2994 c->saved_values = 0; 2995 c->idiag = 0; 2996 c->ssor = 0; 2997 c->keepzeroedrows = a->keepzeroedrows; 2998 c->freedata = PETSC_TRUE; 2999 c->xtoy = 0; 3000 c->XtoY = 0; 3001 3002 c->nz = a->nz; 3003 c->maxnz = a->maxnz; 3004 C->preallocated = PETSC_TRUE; 3005 3006 c->compressedrow.use = a->compressedrow.use; 3007 c->compressedrow.nrows = a->compressedrow.nrows; 3008 c->compressedrow.checked = a->compressedrow.checked; 3009 if ( a->compressedrow.checked && a->compressedrow.use){ 3010 i = a->compressedrow.nrows; 3011 ierr = PetscMalloc((2*i+1)*sizeof(PetscInt),&c->compressedrow.i);CHKERRQ(ierr); 3012 c->compressedrow.rindex = c->compressedrow.i + i + 1; 3013 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3014 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3015 } else { 3016 c->compressedrow.use = PETSC_FALSE; 3017 c->compressedrow.i = PETSC_NULL; 3018 c->compressedrow.rindex = PETSC_NULL; 3019 } 3020 C->same_nonzero = A->same_nonzero; 3021 ierr = MatDuplicate_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3022 3023 *B = C; 3024 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 3025 PetscFunctionReturn(0); 3026 } 3027 3028 #undef __FUNCT__ 3029 #define __FUNCT__ "MatLoad_SeqAIJ" 3030 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, MatType type,Mat *A) 3031 { 3032 Mat_SeqAIJ *a; 3033 Mat B; 3034 PetscErrorCode ierr; 3035 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3036 int fd; 3037 PetscMPIInt size; 3038 MPI_Comm comm; 3039 3040 PetscFunctionBegin; 3041 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3042 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3043 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 3044 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3045 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3046 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3047 M = header[1]; N = header[2]; nz = header[3]; 3048 3049 if (nz < 0) { 3050 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3051 } 3052 3053 /* read in row lengths */ 3054 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3055 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3056 3057 /* check if sum of rowlengths is same as nz */ 3058 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3059 if (sum != nz) SETERRQ2(PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3060 3061 /* create our matrix */ 3062 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3063 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3064 ierr = MatSetType(B,type);CHKERRQ(ierr); 3065 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3066 a = (Mat_SeqAIJ*)B->data; 3067 3068 /* read in column indices and adjust for Fortran indexing*/ 3069 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3070 3071 /* read in nonzero values */ 3072 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3073 3074 /* set matrix "i" values */ 3075 a->i[0] = 0; 3076 for (i=1; i<= M; i++) { 3077 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3078 a->ilen[i-1] = rowlengths[i-1]; 3079 } 3080 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3081 3082 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3083 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3084 *A = B; 3085 PetscFunctionReturn(0); 3086 } 3087 3088 #undef __FUNCT__ 3089 #define __FUNCT__ "MatEqual_SeqAIJ" 3090 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3091 { 3092 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3093 PetscErrorCode ierr; 3094 3095 PetscFunctionBegin; 3096 /* If the matrix dimensions are not equal,or no of nonzeros */ 3097 if ((A->rmap.n != B->rmap.n) || (A->cmap.n != B->cmap.n) ||(a->nz != b->nz)) { 3098 *flg = PETSC_FALSE; 3099 PetscFunctionReturn(0); 3100 } 3101 3102 /* if the a->i are the same */ 3103 ierr = PetscMemcmp(a->i,b->i,(A->rmap.n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3104 if (!*flg) PetscFunctionReturn(0); 3105 3106 /* if a->j are the same */ 3107 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3108 if (!*flg) PetscFunctionReturn(0); 3109 3110 /* if a->a are the same */ 3111 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3112 3113 PetscFunctionReturn(0); 3114 3115 } 3116 3117 #undef __FUNCT__ 3118 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3119 /*@ 3120 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3121 provided by the user. 3122 3123 Collective on MPI_Comm 3124 3125 Input Parameters: 3126 + comm - must be an MPI communicator of size 1 3127 . m - number of rows 3128 . n - number of columns 3129 . i - row indices 3130 . j - column indices 3131 - a - matrix values 3132 3133 Output Parameter: 3134 . mat - the matrix 3135 3136 Level: intermediate 3137 3138 Notes: 3139 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3140 once the matrix is destroyed 3141 3142 You cannot set new nonzero locations into this matrix, that will generate an error. 3143 3144 The i and j indices are 0 based 3145 3146 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3147 3148 @*/ 3149 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3150 { 3151 PetscErrorCode ierr; 3152 PetscInt ii; 3153 Mat_SeqAIJ *aij; 3154 3155 PetscFunctionBegin; 3156 if (i[0]) { 3157 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3158 } 3159 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3160 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3161 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3162 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3163 aij = (Mat_SeqAIJ*)(*mat)->data; 3164 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3165 3166 aij->i = i; 3167 aij->j = j; 3168 aij->a = a; 3169 aij->singlemalloc = PETSC_FALSE; 3170 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3171 aij->freedata = PETSC_FALSE; 3172 3173 for (ii=0; ii<m; ii++) { 3174 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3175 #if defined(PETSC_USE_DEBUG) 3176 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]); 3177 #endif 3178 } 3179 #if defined(PETSC_USE_DEBUG) 3180 for (ii=0; ii<aij->i[m]; ii++) { 3181 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3182 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3183 } 3184 #endif 3185 3186 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3187 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3188 PetscFunctionReturn(0); 3189 } 3190 3191 #undef __FUNCT__ 3192 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3193 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3194 { 3195 PetscErrorCode ierr; 3196 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3197 3198 PetscFunctionBegin; 3199 if (coloring->ctype == IS_COLORING_LOCAL) { 3200 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3201 a->coloring = coloring; 3202 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3203 PetscInt i,*larray; 3204 ISColoring ocoloring; 3205 ISColoringValue *colors; 3206 3207 /* set coloring for diagonal portion */ 3208 ierr = PetscMalloc((A->cmap.n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3209 for (i=0; i<A->cmap.n; i++) { 3210 larray[i] = i; 3211 } 3212 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap.n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3213 ierr = PetscMalloc((A->cmap.n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3214 for (i=0; i<A->cmap.n; i++) { 3215 colors[i] = coloring->colors[larray[i]]; 3216 } 3217 ierr = PetscFree(larray);CHKERRQ(ierr); 3218 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap.n,colors,&ocoloring);CHKERRQ(ierr); 3219 a->coloring = ocoloring; 3220 } 3221 PetscFunctionReturn(0); 3222 } 3223 3224 #if defined(PETSC_HAVE_ADIC) 3225 EXTERN_C_BEGIN 3226 #include "adic/ad_utils.h" 3227 EXTERN_C_END 3228 3229 #undef __FUNCT__ 3230 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3231 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3232 { 3233 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3234 PetscInt m = A->rmap.n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3235 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3236 ISColoringValue *color; 3237 3238 PetscFunctionBegin; 3239 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3240 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3241 color = a->coloring->colors; 3242 /* loop over rows */ 3243 for (i=0; i<m; i++) { 3244 nz = ii[i+1] - ii[i]; 3245 /* loop over columns putting computed value into matrix */ 3246 for (j=0; j<nz; j++) { 3247 *v++ = values[color[*jj++]]; 3248 } 3249 values += nlen; /* jump to next row of derivatives */ 3250 } 3251 PetscFunctionReturn(0); 3252 } 3253 #endif 3254 3255 #undef __FUNCT__ 3256 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3257 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3258 { 3259 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3260 PetscInt m = A->rmap.n,*ii = a->i,*jj = a->j,nz,i,j; 3261 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3262 ISColoringValue *color; 3263 3264 PetscFunctionBegin; 3265 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3266 color = a->coloring->colors; 3267 /* loop over rows */ 3268 for (i=0; i<m; i++) { 3269 nz = ii[i+1] - ii[i]; 3270 /* loop over columns putting computed value into matrix */ 3271 for (j=0; j<nz; j++) { 3272 *v++ = values[color[*jj++]]; 3273 } 3274 values += nl; /* jump to next row of derivatives */ 3275 } 3276 PetscFunctionReturn(0); 3277 } 3278 3279 /* 3280 Special version for direct calls from Fortran 3281 */ 3282 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3283 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3284 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3285 #define matsetvaluesseqaij_ matsetvaluesseqaij 3286 #endif 3287 3288 /* Change these macros so can be used in void function */ 3289 #undef CHKERRQ 3290 #define CHKERRQ(ierr) CHKERRABORT(A->comm,ierr) 3291 #undef SETERRQ2 3292 #define SETERRQ2(ierr,b,c,d) CHKERRABORT(A->comm,ierr) 3293 3294 EXTERN_C_BEGIN 3295 #undef __FUNCT__ 3296 #define __FUNCT__ "matsetvaluesseqaij_" 3297 void PETSCMAT_DLLEXPORT matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis) 3298 { 3299 Mat A = *AA; 3300 PetscInt m = *mm, n = *nn; 3301 InsertMode is = *isis; 3302 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3303 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3304 PetscInt *imax,*ai,*ailen; 3305 PetscErrorCode ierr; 3306 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3307 PetscScalar *ap,value,*aa; 3308 PetscTruth ignorezeroentries = ((a->ignorezeroentries && is == ADD_VALUES) ? PETSC_TRUE:PETSC_FALSE); 3309 PetscTruth roworiented = a->roworiented; 3310 3311 PetscFunctionBegin; 3312 MatPreallocated(A); 3313 imax = a->imax; 3314 ai = a->i; 3315 ailen = a->ilen; 3316 aj = a->j; 3317 aa = a->a; 3318 3319 for (k=0; k<m; k++) { /* loop over added rows */ 3320 row = im[k]; 3321 if (row < 0) continue; 3322 #if defined(PETSC_USE_DEBUG) 3323 if (row >= A->rmap.n) SETERRABORT(A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3324 #endif 3325 rp = aj + ai[row]; ap = aa + ai[row]; 3326 rmax = imax[row]; nrow = ailen[row]; 3327 low = 0; 3328 high = nrow; 3329 for (l=0; l<n; l++) { /* loop over added columns */ 3330 if (in[l] < 0) continue; 3331 #if defined(PETSC_USE_DEBUG) 3332 if (in[l] >= A->cmap.n) SETERRABORT(A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3333 #endif 3334 col = in[l]; 3335 if (roworiented) { 3336 value = v[l + k*n]; 3337 } else { 3338 value = v[k + l*m]; 3339 } 3340 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3341 3342 if (col <= lastcol) low = 0; else high = nrow; 3343 lastcol = col; 3344 while (high-low > 5) { 3345 t = (low+high)/2; 3346 if (rp[t] > col) high = t; 3347 else low = t; 3348 } 3349 for (i=low; i<high; i++) { 3350 if (rp[i] > col) break; 3351 if (rp[i] == col) { 3352 if (is == ADD_VALUES) ap[i] += value; 3353 else ap[i] = value; 3354 goto noinsert; 3355 } 3356 } 3357 if (value == 0.0 && ignorezeroentries) goto noinsert; 3358 if (nonew == 1) goto noinsert; 3359 if (nonew == -1) SETERRABORT(A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3360 MatSeqXAIJReallocateAIJ(a,1,nrow,row,col,rmax,aa,ai,aj,A->rmap.n,rp,ap,imax,nonew); 3361 N = nrow++ - 1; a->nz++; high++; 3362 /* shift up all the later entries in this row */ 3363 for (ii=N; ii>=i; ii--) { 3364 rp[ii+1] = rp[ii]; 3365 ap[ii+1] = ap[ii]; 3366 } 3367 rp[i] = col; 3368 ap[i] = value; 3369 noinsert:; 3370 low = i + 1; 3371 } 3372 ailen[row] = nrow; 3373 } 3374 A->same_nonzero = PETSC_FALSE; 3375 PetscFunctionReturnVoid(); 3376 } 3377 EXTERN_C_END 3378