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