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