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