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