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