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