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