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