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