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