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