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