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