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