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