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