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