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