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