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