1 /*$Id: sbaij.c,v 1.62 2001/08/07 03:03:01 balay Exp $*/ 2 3 /* 4 Defines the basic matrix operations for the BAIJ (compressed row) 5 matrix storage format. 6 */ 7 #include "src/mat/impls/baij/seq/baij.h" /*I "petscmat.h" I*/ 8 #include "src/vec/vecimpl.h" 9 #include "src/inline/spops.h" 10 #include "src/mat/impls/sbaij/seq/sbaij.h" 11 12 #define CHUNKSIZE 10 13 14 /* 15 Checks for missing diagonals 16 */ 17 #undef __FUNCT__ 18 #define __FUNCT__ "MatMissingDiagonal_SeqSBAIJ" 19 int MatMissingDiagonal_SeqSBAIJ(Mat A) 20 { 21 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 22 int *diag,*jj = a->j,i,ierr; 23 24 PetscFunctionBegin; 25 ierr = MatMarkDiagonal_SeqSBAIJ(A);CHKERRQ(ierr); 26 diag = a->diag; 27 for (i=0; i<a->mbs; i++) { 28 if (jj[diag[i]] != i) { 29 SETERRQ1(1,"Matrix is missing diagonal number %d",i); 30 } 31 } 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "MatMarkDiagonal_SeqSBAIJ" 37 int MatMarkDiagonal_SeqSBAIJ(Mat A) 38 { 39 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 40 int i,j,*diag,m = a->mbs,ierr; 41 42 PetscFunctionBegin; 43 if (a->diag) PetscFunctionReturn(0); 44 45 ierr = PetscMalloc((m+1)*sizeof(int),&diag);CHKERRQ(ierr); 46 PetscLogObjectMemory(A,(m+1)*sizeof(int)); 47 for (i=0; i<m; i++) { 48 diag[i] = a->i[i+1]; 49 for (j=a->i[i]; j<a->i[i+1]; j++) { 50 if (a->j[j] == i) { 51 diag[i] = j; 52 break; 53 } 54 } 55 } 56 a->diag = diag; 57 PetscFunctionReturn(0); 58 } 59 60 61 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 62 63 #undef __FUNCT__ 64 #define __FUNCT__ "MatGetRowIJ_SeqSBAIJ" 65 static int MatGetRowIJ_SeqSBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done) 66 { 67 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 68 69 PetscFunctionBegin; 70 if (ia) { 71 SETERRQ(1,"Function not yet written for SBAIJ format, only supports natural ordering"); 72 } 73 *nn = a->mbs; 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "MatRestoreRowIJ_SeqSBAIJ" 79 static int MatRestoreRowIJ_SeqSBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done) 80 { 81 PetscFunctionBegin; 82 if (!ia) PetscFunctionReturn(0); 83 SETERRQ(1,"Function not yet written for SBAIJ format, only supports natural ordering"); 84 /* PetscFunctionReturn(0); */ 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "MatGetBlockSize_SeqSBAIJ" 89 int MatGetBlockSize_SeqSBAIJ(Mat mat,int *bs) 90 { 91 Mat_SeqSBAIJ *sbaij = (Mat_SeqSBAIJ*)mat->data; 92 93 PetscFunctionBegin; 94 *bs = sbaij->bs; 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "MatDestroy_SeqSBAIJ" 100 int MatDestroy_SeqSBAIJ(Mat A) 101 { 102 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 103 int ierr; 104 105 PetscFunctionBegin; 106 #if defined(PETSC_USE_LOG) 107 PetscLogObjectState((PetscObject)A,"Rows=%d, s_NZ=%d",A->m,a->s_nz); 108 #endif 109 ierr = PetscFree(a->a);CHKERRQ(ierr); 110 if (!a->singlemalloc) { 111 ierr = PetscFree(a->i);CHKERRQ(ierr); 112 ierr = PetscFree(a->j);CHKERRQ(ierr); 113 } 114 if (a->row) { 115 ierr = ISDestroy(a->row);CHKERRQ(ierr); 116 } 117 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 118 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 119 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 120 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 121 if (a->mult_work) {ierr = PetscFree(a->mult_work);CHKERRQ(ierr);} 122 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 123 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 124 125 if (a->inew){ 126 ierr = PetscFree(a->inew);CHKERRQ(ierr); 127 a->inew = 0; 128 } 129 ierr = PetscFree(a);CHKERRQ(ierr); 130 PetscFunctionReturn(0); 131 } 132 133 #undef __FUNCT__ 134 #define __FUNCT__ "MatSetOption_SeqSBAIJ" 135 int MatSetOption_SeqSBAIJ(Mat A,MatOption op) 136 { 137 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 138 139 PetscFunctionBegin; 140 switch (op) { 141 case MAT_ROW_ORIENTED: 142 a->roworiented = PETSC_TRUE; 143 break; 144 case MAT_COLUMN_ORIENTED: 145 a->roworiented = PETSC_FALSE; 146 break; 147 case MAT_COLUMNS_SORTED: 148 a->sorted = PETSC_TRUE; 149 break; 150 case MAT_COLUMNS_UNSORTED: 151 a->sorted = PETSC_FALSE; 152 break; 153 case MAT_KEEP_ZEROED_ROWS: 154 a->keepzeroedrows = PETSC_TRUE; 155 break; 156 case MAT_NO_NEW_NONZERO_LOCATIONS: 157 a->nonew = 1; 158 break; 159 case MAT_NEW_NONZERO_LOCATION_ERR: 160 a->nonew = -1; 161 break; 162 case MAT_NEW_NONZERO_ALLOCATION_ERR: 163 a->nonew = -2; 164 break; 165 case MAT_YES_NEW_NONZERO_LOCATIONS: 166 a->nonew = 0; 167 break; 168 case MAT_ROWS_SORTED: 169 case MAT_ROWS_UNSORTED: 170 case MAT_YES_NEW_DIAGONALS: 171 case MAT_IGNORE_OFF_PROC_ENTRIES: 172 case MAT_USE_HASH_TABLE: 173 case MAT_USE_SINGLE_PRECISION_SOLVES: 174 PetscLogInfo(A,"MatSetOption_SeqSBAIJ:Option ignored\n"); 175 break; 176 case MAT_NO_NEW_DIAGONALS: 177 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 178 default: 179 SETERRQ(PETSC_ERR_SUP,"unknown option"); 180 } 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "MatGetRow_SeqSBAIJ" 186 int MatGetRow_SeqSBAIJ(Mat A,int row,int *ncols,int **cols,PetscScalar **v) 187 { 188 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 189 int itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*cols_i,bs2,ierr; 190 MatScalar *aa,*aa_i; 191 PetscScalar *v_i; 192 193 PetscFunctionBegin; 194 bs = a->bs; 195 ai = a->i; 196 aj = a->j; 197 aa = a->a; 198 bs2 = a->bs2; 199 200 if (row < 0 || row >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row out of range"); 201 202 bn = row/bs; /* Block number */ 203 bp = row % bs; /* Block position */ 204 M = ai[bn+1] - ai[bn]; 205 *ncols = bs*M; 206 207 if (v) { 208 *v = 0; 209 if (*ncols) { 210 ierr = PetscMalloc((*ncols+row)*sizeof(PetscScalar),v);CHKERRQ(ierr); 211 for (i=0; i<M; i++) { /* for each block in the block row */ 212 v_i = *v + i*bs; 213 aa_i = aa + bs2*(ai[bn] + i); 214 for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];} 215 } 216 } 217 } 218 219 if (cols) { 220 *cols = 0; 221 if (*ncols) { 222 ierr = PetscMalloc((*ncols+row)*sizeof(int),cols);CHKERRQ(ierr); 223 for (i=0; i<M; i++) { /* for each block in the block row */ 224 cols_i = *cols + i*bs; 225 itmp = bs*aj[ai[bn] + i]; 226 for (j=0; j<bs; j++) {cols_i[j] = itmp++;} 227 } 228 } 229 } 230 231 /*search column A(0:row-1,row) (=A(row,0:row-1)). Could be expensive! */ 232 /* this segment is currently removed, so only entries in the upper triangle are obtained */ 233 #ifdef column_search 234 v_i = *v + M*bs; 235 cols_i = *cols + M*bs; 236 for (i=0; i<bn; i++){ /* for each block row */ 237 M = ai[i+1] - ai[i]; 238 for (j=0; j<M; j++){ 239 itmp = aj[ai[i] + j]; /* block column value */ 240 if (itmp == bn){ 241 aa_i = aa + bs2*(ai[i] + j) + bs*bp; 242 for (k=0; k<bs; k++) { 243 *cols_i++ = i*bs+k; 244 *v_i++ = aa_i[k]; 245 } 246 *ncols += bs; 247 break; 248 } 249 } 250 } 251 #endif 252 253 PetscFunctionReturn(0); 254 } 255 256 #undef __FUNCT__ 257 #define __FUNCT__ "MatRestoreRow_SeqSBAIJ" 258 int MatRestoreRow_SeqSBAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 259 { 260 int ierr; 261 262 PetscFunctionBegin; 263 if (idx) {if (*idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 264 if (v) {if (*v) {ierr = PetscFree(*v);CHKERRQ(ierr);}} 265 PetscFunctionReturn(0); 266 } 267 268 #undef __FUNCT__ 269 #define __FUNCT__ "MatTranspose_SeqSBAIJ" 270 int MatTranspose_SeqSBAIJ(Mat A,Mat *B) 271 { 272 PetscFunctionBegin; 273 SETERRQ(1,"Matrix is symmetric. MatTranspose() should not be called"); 274 /* PetscFunctionReturn(0); */ 275 } 276 277 #undef __FUNCT__ 278 #define __FUNCT__ "MatView_SeqSBAIJ_Binary" 279 static int MatView_SeqSBAIJ_Binary(Mat A,PetscViewer viewer) 280 { 281 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 282 int i,fd,*col_lens,ierr,bs = a->bs,count,*jj,j,k,l,bs2=a->bs2; 283 PetscScalar *aa; 284 FILE *file; 285 286 PetscFunctionBegin; 287 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 288 ierr = PetscMalloc((4+A->m)*sizeof(int),&col_lens);CHKERRQ(ierr); 289 col_lens[0] = MAT_FILE_COOKIE; 290 291 col_lens[1] = A->m; 292 col_lens[2] = A->m; 293 col_lens[3] = a->s_nz*bs2; 294 295 /* store lengths of each row and write (including header) to file */ 296 count = 0; 297 for (i=0; i<a->mbs; i++) { 298 for (j=0; j<bs; j++) { 299 col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]); 300 } 301 } 302 ierr = PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,1);CHKERRQ(ierr); 303 ierr = PetscFree(col_lens);CHKERRQ(ierr); 304 305 /* store column indices (zero start index) */ 306 ierr = PetscMalloc((a->s_nz+1)*bs2*sizeof(int),&jj);CHKERRQ(ierr); 307 count = 0; 308 for (i=0; i<a->mbs; i++) { 309 for (j=0; j<bs; j++) { 310 for (k=a->i[i]; k<a->i[i+1]; k++) { 311 for (l=0; l<bs; l++) { 312 jj[count++] = bs*a->j[k] + l; 313 } 314 } 315 } 316 } 317 ierr = PetscBinaryWrite(fd,jj,bs2*a->s_nz,PETSC_INT,0);CHKERRQ(ierr); 318 ierr = PetscFree(jj);CHKERRQ(ierr); 319 320 /* store nonzero values */ 321 ierr = PetscMalloc((a->s_nz+1)*bs2*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 322 count = 0; 323 for (i=0; i<a->mbs; i++) { 324 for (j=0; j<bs; j++) { 325 for (k=a->i[i]; k<a->i[i+1]; k++) { 326 for (l=0; l<bs; l++) { 327 aa[count++] = a->a[bs2*k + l*bs + j]; 328 } 329 } 330 } 331 } 332 ierr = PetscBinaryWrite(fd,aa,bs2*a->s_nz,PETSC_SCALAR,0);CHKERRQ(ierr); 333 ierr = PetscFree(aa);CHKERRQ(ierr); 334 335 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 336 if (file) { 337 fprintf(file,"-matload_block_size %d\n",a->bs); 338 } 339 PetscFunctionReturn(0); 340 } 341 342 #undef __FUNCT__ 343 #define __FUNCT__ "MatView_SeqSBAIJ_ASCII" 344 static int MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 345 { 346 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 347 int ierr,i,j,bs = a->bs,k,l,bs2=a->bs2; 348 char *name; 349 PetscViewerFormat format; 350 351 PetscFunctionBegin; 352 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 353 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 354 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_LONG) { 355 ierr = PetscViewerASCIIPrintf(viewer," block size is %d\n",bs);CHKERRQ(ierr); 356 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 357 SETERRQ(PETSC_ERR_SUP,"Matlab format not supported"); 358 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 359 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 360 for (i=0; i<a->mbs; i++) { 361 for (j=0; j<bs; j++) { 362 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 363 for (k=a->i[i]; k<a->i[i+1]; k++) { 364 for (l=0; l<bs; l++) { 365 #if defined(PETSC_USE_COMPLEX) 366 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 367 ierr = PetscViewerASCIIPrintf(viewer," %d %g + %g i",bs*a->j[k]+l, 368 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 369 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 370 ierr = PetscViewerASCIIPrintf(viewer," %d %g - %g i",bs*a->j[k]+l, 371 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 372 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 373 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 374 } 375 #else 376 if (a->a[bs2*k + l*bs + j] != 0.0) { 377 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 378 } 379 #endif 380 } 381 } 382 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 383 } 384 } 385 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 386 } else { 387 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 388 for (i=0; i<a->mbs; i++) { 389 for (j=0; j<bs; j++) { 390 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 391 for (k=a->i[i]; k<a->i[i+1]; k++) { 392 for (l=0; l<bs; l++) { 393 #if defined(PETSC_USE_COMPLEX) 394 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 395 ierr = PetscViewerASCIIPrintf(viewer," %d %g + %g i",bs*a->j[k]+l, 396 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 397 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 398 ierr = PetscViewerASCIIPrintf(viewer," %d %g - %g i",bs*a->j[k]+l, 399 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 400 } else { 401 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 402 } 403 #else 404 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 405 #endif 406 } 407 } 408 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 409 } 410 } 411 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 412 } 413 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 414 PetscFunctionReturn(0); 415 } 416 417 #undef __FUNCT__ 418 #define __FUNCT__ "MatView_SeqSBAIJ_Draw_Zoom" 419 static int MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 420 { 421 Mat A = (Mat) Aa; 422 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 423 int row,ierr,i,j,k,l,mbs=a->mbs,color,bs=a->bs,bs2=a->bs2,rank; 424 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 425 MatScalar *aa; 426 MPI_Comm comm; 427 PetscViewer viewer; 428 429 PetscFunctionBegin; 430 /* 431 This is nasty. If this is called from an originally parallel matrix 432 then all processes call this,but only the first has the matrix so the 433 rest should return immediately. 434 */ 435 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 436 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 437 if (rank) PetscFunctionReturn(0); 438 439 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 440 441 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 442 PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric"); 443 444 /* loop over matrix elements drawing boxes */ 445 color = PETSC_DRAW_BLUE; 446 for (i=0,row=0; i<mbs; i++,row+=bs) { 447 for (j=a->i[i]; j<a->i[i+1]; j++) { 448 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 449 x_l = a->j[j]*bs; x_r = x_l + 1.0; 450 aa = a->a + j*bs2; 451 for (k=0; k<bs; k++) { 452 for (l=0; l<bs; l++) { 453 if (PetscRealPart(*aa++) >= 0.) continue; 454 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 455 } 456 } 457 } 458 } 459 color = PETSC_DRAW_CYAN; 460 for (i=0,row=0; i<mbs; i++,row+=bs) { 461 for (j=a->i[i]; j<a->i[i+1]; j++) { 462 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 463 x_l = a->j[j]*bs; x_r = x_l + 1.0; 464 aa = a->a + j*bs2; 465 for (k=0; k<bs; k++) { 466 for (l=0; l<bs; l++) { 467 if (PetscRealPart(*aa++) != 0.) continue; 468 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 469 } 470 } 471 } 472 } 473 474 color = PETSC_DRAW_RED; 475 for (i=0,row=0; i<mbs; i++,row+=bs) { 476 for (j=a->i[i]; j<a->i[i+1]; j++) { 477 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 478 x_l = a->j[j]*bs; x_r = x_l + 1.0; 479 aa = a->a + j*bs2; 480 for (k=0; k<bs; k++) { 481 for (l=0; l<bs; l++) { 482 if (PetscRealPart(*aa++) <= 0.) continue; 483 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 484 } 485 } 486 } 487 } 488 PetscFunctionReturn(0); 489 } 490 491 #undef __FUNCT__ 492 #define __FUNCT__ "MatView_SeqSBAIJ_Draw" 493 static int MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 494 { 495 int ierr; 496 PetscReal xl,yl,xr,yr,w,h; 497 PetscDraw draw; 498 PetscTruth isnull; 499 500 PetscFunctionBegin; 501 502 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 503 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 504 505 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 506 xr = A->m; yr = A->m; h = yr/10.0; w = xr/10.0; 507 xr += w; yr += h; xl = -w; yl = -h; 508 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 509 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 510 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 511 PetscFunctionReturn(0); 512 } 513 514 #undef __FUNCT__ 515 #define __FUNCT__ "MatView_SeqSBAIJ" 516 int MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 517 { 518 int ierr; 519 PetscTruth issocket,isascii,isbinary,isdraw; 520 521 PetscFunctionBegin; 522 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 523 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 524 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 525 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 526 if (issocket) { 527 SETERRQ(PETSC_ERR_SUP,"Socket viewer not supported"); 528 } else if (isascii){ 529 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 530 } else if (isbinary) { 531 ierr = MatView_SeqSBAIJ_Binary(A,viewer);CHKERRQ(ierr); 532 } else if (isdraw) { 533 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 534 } else { 535 SETERRQ1(1,"Viewer type %s not supported by SeqBAIJ matrices",((PetscObject)viewer)->type_name); 536 } 537 PetscFunctionReturn(0); 538 } 539 540 541 #undef __FUNCT__ 542 #define __FUNCT__ "MatGetValues_SeqSBAIJ" 543 int MatGetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *v) 544 { 545 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 546 int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 547 int *ai = a->i,*ailen = a->ilen; 548 int brow,bcol,ridx,cidx,bs=a->bs,bs2=a->bs2; 549 MatScalar *ap,*aa = a->a,zero = 0.0; 550 551 PetscFunctionBegin; 552 for (k=0; k<m; k++) { /* loop over rows */ 553 row = im[k]; brow = row/bs; 554 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative row"); 555 if (row >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 556 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 557 nrow = ailen[brow]; 558 for (l=0; l<n; l++) { /* loop over columns */ 559 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative column"); 560 if (in[l] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 561 col = in[l] ; 562 bcol = col/bs; 563 cidx = col%bs; 564 ridx = row%bs; 565 high = nrow; 566 low = 0; /* assume unsorted */ 567 while (high-low > 5) { 568 t = (low+high)/2; 569 if (rp[t] > bcol) high = t; 570 else low = t; 571 } 572 for (i=low; i<high; i++) { 573 if (rp[i] > bcol) break; 574 if (rp[i] == bcol) { 575 *v++ = ap[bs2*i+bs*cidx+ridx]; 576 goto finished; 577 } 578 } 579 *v++ = zero; 580 finished:; 581 } 582 } 583 PetscFunctionReturn(0); 584 } 585 586 587 #undef __FUNCT__ 588 #define __FUNCT__ "MatSetValuesBlocked_SeqSBAIJ" 589 int MatSetValuesBlocked_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *v,InsertMode is) 590 { 591 PetscFunctionBegin; 592 SETERRQ(1,"Function not yet written for SBAIJ format"); 593 /* PetscFunctionReturn(0); */ 594 } 595 596 #undef __FUNCT__ 597 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 598 int MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 599 { 600 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 601 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 602 int m = A->m,*ip,N,*ailen = a->ilen; 603 int mbs = a->mbs,bs2 = a->bs2,rmax = 0,ierr; 604 MatScalar *aa = a->a,*ap; 605 606 PetscFunctionBegin; 607 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 608 609 if (m) rmax = ailen[0]; 610 for (i=1; i<mbs; i++) { 611 /* move each row back by the amount of empty slots (fshift) before it*/ 612 fshift += imax[i-1] - ailen[i-1]; 613 rmax = PetscMax(rmax,ailen[i]); 614 if (fshift) { 615 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 616 N = ailen[i]; 617 for (j=0; j<N; j++) { 618 ip[j-fshift] = ip[j]; 619 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 620 } 621 } 622 ai[i] = ai[i-1] + ailen[i-1]; 623 } 624 if (mbs) { 625 fshift += imax[mbs-1] - ailen[mbs-1]; 626 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 627 } 628 /* reset ilen and imax for each row */ 629 for (i=0; i<mbs; i++) { 630 ailen[i] = imax[i] = ai[i+1] - ai[i]; 631 } 632 a->s_nz = ai[mbs]; 633 634 /* diagonals may have moved, so kill the diagonal pointers */ 635 if (fshift && a->diag) { 636 ierr = PetscFree(a->diag);CHKERRQ(ierr); 637 PetscLogObjectMemory(A,-(m+1)*sizeof(int)); 638 a->diag = 0; 639 } 640 PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Matrix size: %d X %d, block size %d; storage space: %d unneeded, %d used\n", 641 m,A->m,a->bs,fshift*bs2,a->s_nz*bs2); 642 PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Number of mallocs during MatSetValues is %d\n", 643 a->reallocs); 644 PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Most nonzeros blocks in any row is %d\n",rmax); 645 a->reallocs = 0; 646 A->info.nz_unneeded = (PetscReal)fshift*bs2; 647 648 PetscFunctionReturn(0); 649 } 650 651 /* 652 This function returns an array of flags which indicate the locations of contiguous 653 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 654 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 655 Assume: sizes should be long enough to hold all the values. 656 */ 657 #undef __FUNCT__ 658 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 659 static int MatZeroRows_SeqSBAIJ_Check_Blocks(int idx[],int n,int bs,int sizes[], int *bs_max) 660 { 661 int i,j,k,row; 662 PetscTruth flg; 663 664 PetscFunctionBegin; 665 for (i=0,j=0; i<n; j++) { 666 row = idx[i]; 667 if (row%bs!=0) { /* Not the begining of a block */ 668 sizes[j] = 1; 669 i++; 670 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 671 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 672 i++; 673 } else { /* Begining of the block, so check if the complete block exists */ 674 flg = PETSC_TRUE; 675 for (k=1; k<bs; k++) { 676 if (row+k != idx[i+k]) { /* break in the block */ 677 flg = PETSC_FALSE; 678 break; 679 } 680 } 681 if (flg == PETSC_TRUE) { /* No break in the bs */ 682 sizes[j] = bs; 683 i+= bs; 684 } else { 685 sizes[j] = 1; 686 i++; 687 } 688 } 689 } 690 *bs_max = j; 691 PetscFunctionReturn(0); 692 } 693 694 #undef __FUNCT__ 695 #define __FUNCT__ "MatZeroRows_SeqSBAIJ" 696 int MatZeroRows_SeqSBAIJ(Mat A,IS is,PetscScalar *diag) 697 { 698 Mat_SeqSBAIJ *sbaij=(Mat_SeqSBAIJ*)A->data; 699 int ierr,i,j,k,count,is_n,*is_idx,*rows; 700 int bs=sbaij->bs,bs2=sbaij->bs2,*sizes,row,bs_max; 701 PetscScalar zero = 0.0; 702 MatScalar *aa; 703 704 PetscFunctionBegin; 705 /* Make a copy of the IS and sort it */ 706 ierr = ISGetSize(is,&is_n);CHKERRQ(ierr); 707 ierr = ISGetIndices(is,&is_idx);CHKERRQ(ierr); 708 709 /* allocate memory for rows,sizes */ 710 ierr = PetscMalloc((3*is_n+1)*sizeof(int),&rows);CHKERRQ(ierr); 711 sizes = rows + is_n; 712 713 /* initialize copy IS values to rows, and sort them */ 714 for (i=0; i<is_n; i++) { rows[i] = is_idx[i]; } 715 ierr = PetscSortInt(is_n,rows);CHKERRQ(ierr); 716 if (sbaij->keepzeroedrows) { /* do not change nonzero structure */ 717 for (i=0; i<is_n; i++) { sizes[i] = 1; } /* sizes: size of blocks, = 1 or bs */ 718 bs_max = is_n; /* bs_max: num. of contiguous block row in the row */ 719 } else { 720 ierr = MatZeroRows_SeqSBAIJ_Check_Blocks(rows,is_n,bs,sizes,&bs_max);CHKERRQ(ierr); 721 } 722 ierr = ISRestoreIndices(is,&is_idx);CHKERRQ(ierr); 723 724 for (i=0,j=0; i<bs_max; j+=sizes[i],i++) { 725 row = rows[j]; /* row to be zeroed */ 726 if (row < 0 || row > A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %d out of range",row); 727 count = (sbaij->i[row/bs +1] - sbaij->i[row/bs])*bs; /* num. of elements in the row */ 728 aa = sbaij->a + sbaij->i[row/bs]*bs2 + (row%bs); 729 if (sizes[i] == bs && !sbaij->keepzeroedrows) { 730 if (diag) { 731 if (sbaij->ilen[row/bs] > 0) { 732 sbaij->ilen[row/bs] = 1; 733 sbaij->j[sbaij->i[row/bs]] = row/bs; 734 ierr = PetscMemzero(aa,count*bs*sizeof(MatScalar));CHKERRQ(ierr); 735 } 736 /* Now insert all the diagoanl values for this bs */ 737 for (k=0; k<bs; k++) { 738 ierr = (*A->ops->setvalues)(A,1,rows+j+k,1,rows+j+k,diag,INSERT_VALUES);CHKERRQ(ierr); 739 } 740 } else { /* (!diag) */ 741 sbaij->ilen[row/bs] = 0; 742 } /* end (!diag) */ 743 } else { /* (sizes[i] != bs), broken block */ 744 #if defined (PETSC_USE_DEBUG) 745 if (sizes[i] != 1) SETERRQ(1,"Internal Error. Value should be 1"); 746 #endif 747 for (k=0; k<count; k++) { 748 aa[0] = zero; 749 aa+=bs; 750 } 751 if (diag) { 752 ierr = (*A->ops->setvalues)(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); 753 } 754 } 755 } 756 757 ierr = PetscFree(rows);CHKERRQ(ierr); 758 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 759 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 760 PetscFunctionReturn(0); 761 } 762 763 /* Only add/insert a(i,j) with i<=j (blocks). 764 Any a(i,j) with i>j input by user is ingored. 765 */ 766 767 #undef __FUNCT__ 768 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 769 int MatSetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *v,InsertMode is) 770 { 771 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 772 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 773 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 774 int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol; 775 int ridx,cidx,bs2=a->bs2,ierr; 776 MatScalar *ap,value,*aa=a->a,*bap; 777 778 PetscFunctionBegin; 779 780 for (k=0; k<m; k++) { /* loop over added rows */ 781 row = im[k]; /* row number */ 782 brow = row/bs; /* block row number */ 783 if (row < 0) continue; 784 #if defined(PETSC_USE_BOPT_g) 785 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",row,A->m); 786 #endif 787 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 788 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 789 rmax = imax[brow]; /* maximum space allocated for this row */ 790 nrow = ailen[brow]; /* actual length of this row */ 791 low = 0; 792 793 for (l=0; l<n; l++) { /* loop over added columns */ 794 if (in[l] < 0) continue; 795 #if defined(PETSC_USE_BOPT_g) 796 if (in[l] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",in[l],A->m); 797 #endif 798 col = in[l]; 799 bcol = col/bs; /* block col number */ 800 801 if (brow <= bcol){ 802 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 803 if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ 804 /* element value a(k,l) */ 805 if (roworiented) { 806 value = v[l + k*n]; 807 } else { 808 value = v[k + l*m]; 809 } 810 811 /* move pointer bap to a(k,l) quickly and add/insert value */ 812 if (!sorted) low = 0; high = nrow; 813 while (high-low > 7) { 814 t = (low+high)/2; 815 if (rp[t] > bcol) high = t; 816 else low = t; 817 } 818 for (i=low; i<high; i++) { 819 /* printf("The loop of i=low.., rp[%d]=%d\n",i,rp[i]); */ 820 if (rp[i] > bcol) break; 821 if (rp[i] == bcol) { 822 bap = ap + bs2*i + bs*cidx + ridx; 823 if (is == ADD_VALUES) *bap += value; 824 else *bap = value; 825 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 826 if (brow == bcol && ridx < cidx){ 827 bap = ap + bs2*i + bs*ridx + cidx; 828 if (is == ADD_VALUES) *bap += value; 829 else *bap = value; 830 } 831 goto noinsert1; 832 } 833 } 834 835 if (nonew == 1) goto noinsert1; 836 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 837 if (nrow >= rmax) { 838 /* there is no extra room in row, therefore enlarge */ 839 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 840 MatScalar *new_a; 841 842 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 843 844 /* Malloc new storage space */ 845 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 846 ierr = PetscMalloc(len,&new_a);CHKERRQ(ierr); 847 new_j = (int*)(new_a + bs2*new_nz); 848 new_i = new_j + new_nz; 849 850 /* copy over old data into new slots */ 851 for (ii=0; ii<brow+1; ii++) {new_i[ii] = ai[ii];} 852 for (ii=brow+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 853 ierr = PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int));CHKERRQ(ierr); 854 len = (new_nz - CHUNKSIZE - ai[brow] - nrow); 855 ierr = PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow,len*sizeof(int));CHKERRQ(ierr); 856 ierr = PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(MatScalar));CHKERRQ(ierr); 857 ierr = PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));CHKERRQ(ierr); 858 ierr = PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE),aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(MatScalar));CHKERRQ(ierr); 859 /* free up old matrix storage */ 860 ierr = PetscFree(a->a);CHKERRQ(ierr); 861 if (!a->singlemalloc) { 862 ierr = PetscFree(a->i);CHKERRQ(ierr); 863 ierr = PetscFree(a->j);CHKERRQ(ierr); 864 } 865 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 866 a->singlemalloc = PETSC_TRUE; 867 868 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 869 rmax = imax[brow] = imax[brow] + CHUNKSIZE; 870 PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 871 a->s_maxnz += bs2*CHUNKSIZE; 872 a->reallocs++; 873 a->s_nz++; 874 } 875 876 N = nrow++ - 1; 877 /* shift up all the later entries in this row */ 878 for (ii=N; ii>=i; ii--) { 879 rp[ii+1] = rp[ii]; 880 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 881 } 882 if (N>=i) { 883 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 884 } 885 rp[i] = bcol; 886 ap[bs2*i + bs*cidx + ridx] = value; 887 noinsert1:; 888 low = i; 889 /* } */ 890 } 891 } /* end of if .. if.. */ 892 } /* end of loop over added columns */ 893 ailen[brow] = nrow; 894 } /* end of loop over added rows */ 895 896 PetscFunctionReturn(0); 897 } 898 899 extern int MatCholeskyFactorSymbolic_SeqSBAIJ(Mat,IS,PetscReal,Mat*); 900 extern int MatCholeskyFactor_SeqSBAIJ(Mat,IS,PetscReal); 901 extern int MatIncreaseOverlap_SeqSBAIJ(Mat,int,IS*,int); 902 extern int MatGetSubMatrix_SeqSBAIJ(Mat,IS,IS,int,MatReuse,Mat*); 903 extern int MatGetSubMatrices_SeqSBAIJ(Mat,int,IS*,IS*,MatReuse,Mat**); 904 extern int MatMultTranspose_SeqSBAIJ(Mat,Vec,Vec); 905 extern int MatMultTransposeAdd_SeqSBAIJ(Mat,Vec,Vec,Vec); 906 extern int MatScale_SeqSBAIJ(PetscScalar*,Mat); 907 extern int MatNorm_SeqSBAIJ(Mat,NormType,PetscReal *); 908 extern int MatEqual_SeqSBAIJ(Mat,Mat,PetscTruth*); 909 extern int MatGetDiagonal_SeqSBAIJ(Mat,Vec); 910 extern int MatDiagonalScale_SeqSBAIJ(Mat,Vec,Vec); 911 extern int MatGetInfo_SeqSBAIJ(Mat,MatInfoType,MatInfo *); 912 extern int MatZeroEntries_SeqSBAIJ(Mat); 913 extern int MatGetRowMax_SeqSBAIJ(Mat,Vec); 914 915 extern int MatSolve_SeqSBAIJ_N(Mat,Vec,Vec); 916 extern int MatSolve_SeqSBAIJ_1(Mat,Vec,Vec); 917 extern int MatSolve_SeqSBAIJ_2(Mat,Vec,Vec); 918 extern int MatSolve_SeqSBAIJ_3(Mat,Vec,Vec); 919 extern int MatSolve_SeqSBAIJ_4(Mat,Vec,Vec); 920 extern int MatSolve_SeqSBAIJ_5(Mat,Vec,Vec); 921 extern int MatSolve_SeqSBAIJ_6(Mat,Vec,Vec); 922 extern int MatSolve_SeqSBAIJ_7(Mat,Vec,Vec); 923 extern int MatSolveTranspose_SeqSBAIJ_7(Mat,Vec,Vec); 924 extern int MatSolveTranspose_SeqSBAIJ_6(Mat,Vec,Vec); 925 extern int MatSolveTranspose_SeqSBAIJ_5(Mat,Vec,Vec); 926 extern int MatSolveTranspose_SeqSBAIJ_4(Mat,Vec,Vec); 927 extern int MatSolveTranspose_SeqSBAIJ_3(Mat,Vec,Vec); 928 extern int MatSolveTranspose_SeqSBAIJ_2(Mat,Vec,Vec); 929 extern int MatSolveTranspose_SeqSBAIJ_1(Mat,Vec,Vec); 930 931 extern int MatSolve_SeqSBAIJ_1_NaturalOrdering(Mat,Vec,Vec); 932 extern int MatSolve_SeqSBAIJ_2_NaturalOrdering(Mat,Vec,Vec); 933 extern int MatSolve_SeqSBAIJ_3_NaturalOrdering(Mat,Vec,Vec); 934 extern int MatSolve_SeqSBAIJ_4_NaturalOrdering(Mat,Vec,Vec); 935 extern int MatSolve_SeqSBAIJ_5_NaturalOrdering(Mat,Vec,Vec); 936 extern int MatSolve_SeqSBAIJ_6_NaturalOrdering(Mat,Vec,Vec); 937 extern int MatSolve_SeqSBAIJ_7_NaturalOrdering(Mat,Vec,Vec); 938 extern int MatSolve_SeqSBAIJ_N_NaturalOrdering(Mat,Vec,Vec); 939 940 extern int MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat,Mat*); 941 extern int MatCholeskyFactorNumeric_SeqSBAIJ_1(Mat,Mat*); 942 extern int MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat,Mat*); 943 extern int MatCholeskyFactorNumeric_SeqSBAIJ_3(Mat,Mat*); 944 extern int MatCholeskyFactorNumeric_SeqSBAIJ_4(Mat,Mat*); 945 extern int MatCholeskyFactorNumeric_SeqSBAIJ_5(Mat,Mat*); 946 extern int MatCholeskyFactorNumeric_SeqSBAIJ_6(Mat,Mat*); 947 extern int MatCholeskyFactorNumeric_SeqSBAIJ_7(Mat,Mat*); 948 949 extern int MatMult_SeqSBAIJ_1(Mat,Vec,Vec); 950 extern int MatMult_SeqSBAIJ_2(Mat,Vec,Vec); 951 extern int MatMult_SeqSBAIJ_3(Mat,Vec,Vec); 952 extern int MatMult_SeqSBAIJ_4(Mat,Vec,Vec); 953 extern int MatMult_SeqSBAIJ_5(Mat,Vec,Vec); 954 extern int MatMult_SeqSBAIJ_6(Mat,Vec,Vec); 955 extern int MatMult_SeqSBAIJ_7(Mat,Vec,Vec); 956 extern int MatMult_SeqSBAIJ_N(Mat,Vec,Vec); 957 958 extern int MatMultAdd_SeqSBAIJ_1(Mat,Vec,Vec,Vec); 959 extern int MatMultAdd_SeqSBAIJ_2(Mat,Vec,Vec,Vec); 960 extern int MatMultAdd_SeqSBAIJ_3(Mat,Vec,Vec,Vec); 961 extern int MatMultAdd_SeqSBAIJ_4(Mat,Vec,Vec,Vec); 962 extern int MatMultAdd_SeqSBAIJ_5(Mat,Vec,Vec,Vec); 963 extern int MatMultAdd_SeqSBAIJ_6(Mat,Vec,Vec,Vec); 964 extern int MatMultAdd_SeqSBAIJ_7(Mat,Vec,Vec,Vec); 965 extern int MatMultAdd_SeqSBAIJ_N(Mat,Vec,Vec,Vec); 966 967 #ifdef HAVE_MatICCFactor 968 /* modefied from MatILUFactor_SeqSBAIJ, needs further work! */ 969 #undef __FUNCT__ 970 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 971 int MatICCFactor_SeqSBAIJ(Mat inA,IS row,PetscReal fill,int level) 972 { 973 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 974 Mat outA; 975 int ierr; 976 PetscTruth row_identity,col_identity; 977 978 PetscFunctionBegin; 979 /* 980 if (level != 0) SETERRQ(PETSC_ERR_SUP,"Only levels = 0 supported for in-place ILU"); 981 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 982 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 983 if (!row_identity || !col_identity) { 984 SETERRQ(1,"Row and column permutations must be identity for in-place ICC"); 985 } 986 */ 987 988 outA = inA; 989 inA->factor = FACTOR_CHOLESKY; 990 991 if (!a->diag) { 992 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 993 } 994 /* 995 Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver 996 for ILU(0) factorization with natural ordering 997 */ 998 switch (a->bs) { 999 case 1: 1000 inA->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering; 1001 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering; 1002 PetscLoginfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering solvetrans BS=1\n"); 1003 case 2: 1004 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1005 inA->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering; 1006 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering; 1007 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=2\n"); 1008 break; 1009 case 3: 1010 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1011 inA->ops->solve = MatSolve_SeqSBAIJ_3_NaturalOrdering; 1012 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_3_NaturalOrdering; 1013 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=3\n"); 1014 break; 1015 case 4: 1016 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1017 inA->ops->solve = MatSolve_SeqSBAIJ_4_NaturalOrdering; 1018 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_4_NaturalOrdering; 1019 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=4\n"); 1020 break; 1021 case 5: 1022 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1023 inA->ops->solve = MatSolve_SeqSBAIJ_5_NaturalOrdering; 1024 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_5_NaturalOrdering; 1025 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=5\n"); 1026 break; 1027 case 6: 1028 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1029 inA->ops->solve = MatSolve_SeqSBAIJ_6_NaturalOrdering; 1030 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_6_NaturalOrdering; 1031 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=6\n"); 1032 break; 1033 case 7: 1034 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1035 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_7_NaturalOrdering; 1036 inA->ops->solve = MatSolve_SeqSBAIJ_7_NaturalOrdering; 1037 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=7\n"); 1038 break; 1039 default: 1040 a->row = row; 1041 a->icol = col; 1042 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1043 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1044 1045 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1046 ierr = ISInvertPermutation(col,PETSC_DECIDE, &(a->icol));CHKERRQ(ierr); 1047 PetscLogObjectParent(inA,a->icol); 1048 1049 if (!a->solve_work) { 1050 ierr = PetscMalloc((A->m+a->bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1051 PetscLogObjectMemory(inA,(A->m+a->bs)*sizeof(PetscScalar)); 1052 } 1053 } 1054 1055 ierr = MatCholeskyFactorNumeric(inA,&outA);CHKERRQ(ierr); 1056 1057 PetscFunctionReturn(0); 1058 } 1059 #endif 1060 1061 #undef __FUNCT__ 1062 #define __FUNCT__ "MatPrintHelp_SeqSBAIJ" 1063 int MatPrintHelp_SeqSBAIJ(Mat A) 1064 { 1065 static PetscTruth called = PETSC_FALSE; 1066 MPI_Comm comm = A->comm; 1067 int ierr; 1068 1069 PetscFunctionBegin; 1070 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1071 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQSBAIJ and MATMPISBAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1072 ierr = (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n");CHKERRQ(ierr); 1073 PetscFunctionReturn(0); 1074 } 1075 1076 EXTERN_C_BEGIN 1077 #undef __FUNCT__ 1078 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1079 int MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,int *indices) 1080 { 1081 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1082 int i,nz,n; 1083 1084 PetscFunctionBegin; 1085 nz = baij->s_maxnz; 1086 n = mat->n; 1087 for (i=0; i<nz; i++) { 1088 baij->j[i] = indices[i]; 1089 } 1090 baij->s_nz = nz; 1091 for (i=0; i<n; i++) { 1092 baij->ilen[i] = baij->imax[i]; 1093 } 1094 1095 PetscFunctionReturn(0); 1096 } 1097 EXTERN_C_END 1098 1099 #undef __FUNCT__ 1100 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1101 /*@ 1102 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1103 in the matrix. 1104 1105 Input Parameters: 1106 + mat - the SeqSBAIJ matrix 1107 - indices - the column indices 1108 1109 Level: advanced 1110 1111 Notes: 1112 This can be called if you have precomputed the nonzero structure of the 1113 matrix and want to provide it to the matrix object to improve the performance 1114 of the MatSetValues() operation. 1115 1116 You MUST have set the correct numbers of nonzeros per row in the call to 1117 MatCreateSeqSBAIJ(). 1118 1119 MUST be called before any calls to MatSetValues() 1120 1121 .seealso: MatCreateSeqSBAIJ 1122 @*/ 1123 int MatSeqSBAIJSetColumnIndices(Mat mat,int *indices) 1124 { 1125 int ierr,(*f)(Mat,int *); 1126 1127 PetscFunctionBegin; 1128 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1129 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 1130 if (f) { 1131 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1132 } else { 1133 SETERRQ(1,"Wrong type of matrix to set column indices"); 1134 } 1135 PetscFunctionReturn(0); 1136 } 1137 1138 #undef __FUNCT__ 1139 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1140 int MatSetUpPreallocation_SeqSBAIJ(Mat A) 1141 { 1142 int ierr; 1143 1144 PetscFunctionBegin; 1145 ierr = MatSeqSBAIJSetPreallocation(A,1,PETSC_DEFAULT,0);CHKERRQ(ierr); 1146 PetscFunctionReturn(0); 1147 } 1148 1149 /* -------------------------------------------------------------------*/ 1150 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1151 MatGetRow_SeqSBAIJ, 1152 MatRestoreRow_SeqSBAIJ, 1153 MatMult_SeqSBAIJ_N, 1154 MatMultAdd_SeqSBAIJ_N, 1155 MatMultTranspose_SeqSBAIJ, 1156 MatMultTransposeAdd_SeqSBAIJ, 1157 MatSolve_SeqSBAIJ_N, 1158 0, 1159 0, 1160 0, 1161 0, 1162 MatCholeskyFactor_SeqSBAIJ, 1163 MatRelax_SeqSBAIJ, 1164 MatTranspose_SeqSBAIJ, 1165 MatGetInfo_SeqSBAIJ, 1166 MatEqual_SeqSBAIJ, 1167 MatGetDiagonal_SeqSBAIJ, 1168 MatDiagonalScale_SeqSBAIJ, 1169 MatNorm_SeqSBAIJ, 1170 0, 1171 MatAssemblyEnd_SeqSBAIJ, 1172 0, 1173 MatSetOption_SeqSBAIJ, 1174 MatZeroEntries_SeqSBAIJ, 1175 MatZeroRows_SeqSBAIJ, 1176 0, 1177 0, 1178 MatCholeskyFactorSymbolic_SeqSBAIJ, 1179 MatCholeskyFactorNumeric_SeqSBAIJ_N, 1180 MatSetUpPreallocation_SeqSBAIJ, 1181 0, 1182 MatICCFactorSymbolic_SeqSBAIJ, 1183 0, 1184 0, 1185 MatDuplicate_SeqSBAIJ, 1186 0, 1187 0, 1188 0, 1189 0, 1190 0, 1191 MatGetSubMatrices_SeqSBAIJ, 1192 MatIncreaseOverlap_SeqSBAIJ, 1193 MatGetValues_SeqSBAIJ, 1194 0, 1195 MatPrintHelp_SeqSBAIJ, 1196 MatScale_SeqSBAIJ, 1197 0, 1198 0, 1199 0, 1200 MatGetBlockSize_SeqSBAIJ, 1201 MatGetRowIJ_SeqSBAIJ, 1202 MatRestoreRowIJ_SeqSBAIJ, 1203 0, 1204 0, 1205 0, 1206 0, 1207 0, 1208 0, 1209 MatSetValuesBlocked_SeqSBAIJ, 1210 MatGetSubMatrix_SeqSBAIJ, 1211 0, 1212 0, 1213 MatGetPetscMaps_Petsc, 1214 0, 1215 0, 1216 0, 1217 0, 1218 0, 1219 0, 1220 MatGetRowMax_SeqSBAIJ}; 1221 1222 EXTERN_C_BEGIN 1223 #undef __FUNCT__ 1224 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1225 int MatStoreValues_SeqSBAIJ(Mat mat) 1226 { 1227 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1228 int nz = aij->i[mat->m]*aij->bs*aij->bs2; 1229 int ierr; 1230 1231 PetscFunctionBegin; 1232 if (aij->nonew != 1) { 1233 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1234 } 1235 1236 /* allocate space for values if not already there */ 1237 if (!aij->saved_values) { 1238 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1239 } 1240 1241 /* copy values over */ 1242 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1243 PetscFunctionReturn(0); 1244 } 1245 EXTERN_C_END 1246 1247 EXTERN_C_BEGIN 1248 #undef __FUNCT__ 1249 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1250 int MatRetrieveValues_SeqSBAIJ(Mat mat) 1251 { 1252 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1253 int nz = aij->i[mat->m]*aij->bs*aij->bs2,ierr; 1254 1255 PetscFunctionBegin; 1256 if (aij->nonew != 1) { 1257 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1258 } 1259 if (!aij->saved_values) { 1260 SETERRQ(1,"Must call MatStoreValues(A);first"); 1261 } 1262 1263 /* copy values over */ 1264 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1265 PetscFunctionReturn(0); 1266 } 1267 EXTERN_C_END 1268 1269 EXTERN_C_BEGIN 1270 #undef __FUNCT__ 1271 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1272 int MatCreate_SeqSBAIJ(Mat B) 1273 { 1274 Mat_SeqSBAIJ *b; 1275 int ierr,size; 1276 1277 PetscFunctionBegin; 1278 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1279 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1280 B->m = B->M = PetscMax(B->m,B->M); 1281 B->n = B->N = PetscMax(B->n,B->N); 1282 1283 ierr = PetscNew(Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1284 B->data = (void*)b; 1285 ierr = PetscMemzero(b,sizeof(Mat_SeqSBAIJ));CHKERRQ(ierr); 1286 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1287 B->ops->destroy = MatDestroy_SeqSBAIJ; 1288 B->ops->view = MatView_SeqSBAIJ; 1289 B->factor = 0; 1290 B->lupivotthreshold = 1.0; 1291 B->mapping = 0; 1292 b->row = 0; 1293 b->icol = 0; 1294 b->reallocs = 0; 1295 b->saved_values = 0; 1296 1297 ierr = PetscMapCreateMPI(B->comm,B->m,B->M,&B->rmap);CHKERRQ(ierr); 1298 ierr = PetscMapCreateMPI(B->comm,B->n,B->N,&B->cmap);CHKERRQ(ierr); 1299 1300 b->sorted = PETSC_FALSE; 1301 b->roworiented = PETSC_TRUE; 1302 b->nonew = 0; 1303 b->diag = 0; 1304 b->solve_work = 0; 1305 b->mult_work = 0; 1306 b->spptr = 0; 1307 b->keepzeroedrows = PETSC_FALSE; 1308 1309 b->inew = 0; 1310 b->jnew = 0; 1311 b->anew = 0; 1312 b->a2anew = 0; 1313 b->permute = PETSC_FALSE; 1314 1315 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1316 "MatStoreValues_SeqSBAIJ", 1317 (void*)MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1318 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1319 "MatRetrieveValues_SeqSBAIJ", 1320 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1321 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1322 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1323 (void*)MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1324 PetscFunctionReturn(0); 1325 } 1326 EXTERN_C_END 1327 1328 #undef __FUNCT__ 1329 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1330 /*@C 1331 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1332 compressed row) format. For good matrix assembly performance the 1333 user should preallocate the matrix storage by setting the parameter nz 1334 (or the array nnz). By setting these parameters accurately, performance 1335 during matrix assembly can be increased by more than a factor of 50. 1336 1337 Collective on Mat 1338 1339 Input Parameters: 1340 + A - the symmetric matrix 1341 . bs - size of block 1342 . nz - number of block nonzeros per block row (same for all rows) 1343 - nnz - array containing the number of block nonzeros in the various block rows 1344 (possibly different for each block row) or PETSC_NULL 1345 1346 Options Database Keys: 1347 . -mat_no_unroll - uses code that does not unroll the loops in the 1348 block calculations (much slower) 1349 . -mat_block_size - size of the blocks to use 1350 1351 Level: intermediate 1352 1353 Notes: 1354 The block AIJ format is compatible with standard Fortran 77 1355 storage. That is, the stored row and column indices can begin at 1356 either one (as in Fortran) or zero. See the users' manual for details. 1357 1358 Specify the preallocated storage with either nz or nnz (not both). 1359 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1360 allocation. For additional details, see the users manual chapter on 1361 matrices. 1362 1363 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1364 @*/ 1365 int MatSeqSBAIJSetPreallocation(Mat B,int bs,int nz,int *nnz) 1366 { 1367 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1368 int i,len,ierr,mbs,bs2; 1369 PetscTruth flg; 1370 int s_nz; 1371 1372 PetscFunctionBegin; 1373 B->preallocated = PETSC_TRUE; 1374 ierr = PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1375 mbs = B->m/bs; 1376 bs2 = bs*bs; 1377 1378 if (mbs*bs != B->m) { 1379 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1380 } 1381 1382 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1383 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 1384 if (nnz) { 1385 for (i=0; i<mbs; i++) { 1386 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 1387 if (nnz[i] > mbs) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %d value %d rowlength %d",i,nnz[i],mbs); 1388 } 1389 } 1390 1391 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1392 if (!flg) { 1393 switch (bs) { 1394 case 1: 1395 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1396 B->ops->solve = MatSolve_SeqSBAIJ_1; 1397 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_1; 1398 B->ops->mult = MatMult_SeqSBAIJ_1; 1399 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1400 break; 1401 case 2: 1402 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1403 B->ops->solve = MatSolve_SeqSBAIJ_2; 1404 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_2; 1405 B->ops->mult = MatMult_SeqSBAIJ_2; 1406 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1407 break; 1408 case 3: 1409 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1410 B->ops->solve = MatSolve_SeqSBAIJ_3; 1411 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_3; 1412 B->ops->mult = MatMult_SeqSBAIJ_3; 1413 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1414 break; 1415 case 4: 1416 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1417 B->ops->solve = MatSolve_SeqSBAIJ_4; 1418 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_4; 1419 B->ops->mult = MatMult_SeqSBAIJ_4; 1420 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1421 break; 1422 case 5: 1423 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1424 B->ops->solve = MatSolve_SeqSBAIJ_5; 1425 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_5; 1426 B->ops->mult = MatMult_SeqSBAIJ_5; 1427 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1428 break; 1429 case 6: 1430 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1431 B->ops->solve = MatSolve_SeqSBAIJ_6; 1432 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_6; 1433 B->ops->mult = MatMult_SeqSBAIJ_6; 1434 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1435 break; 1436 case 7: 1437 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1438 B->ops->solve = MatSolve_SeqSBAIJ_7; 1439 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_7; 1440 B->ops->mult = MatMult_SeqSBAIJ_7; 1441 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1442 break; 1443 } 1444 } 1445 1446 b->mbs = mbs; 1447 b->nbs = mbs; 1448 ierr = PetscMalloc((mbs+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 1449 if (!nnz) { 1450 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1451 else if (nz <= 0) nz = 1; 1452 for (i=0; i<mbs; i++) { 1453 b->imax[i] = nz; 1454 } 1455 nz = nz*mbs; /* total nz */ 1456 } else { 1457 nz = 0; 1458 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1459 } 1460 /* s_nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1461 s_nz = nz; 1462 1463 /* allocate the matrix space */ 1464 len = s_nz*sizeof(int) + s_nz*bs2*sizeof(MatScalar) + (B->m+1)*sizeof(int); 1465 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 1466 ierr = PetscMemzero(b->a,s_nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1467 b->j = (int*)(b->a + s_nz*bs2); 1468 ierr = PetscMemzero(b->j,s_nz*sizeof(int));CHKERRQ(ierr); 1469 b->i = b->j + s_nz; 1470 b->singlemalloc = PETSC_TRUE; 1471 1472 /* pointer to beginning of each row */ 1473 b->i[0] = 0; 1474 for (i=1; i<mbs+1; i++) { 1475 b->i[i] = b->i[i-1] + b->imax[i-1]; 1476 } 1477 1478 /* b->ilen will count nonzeros in each block row so far. */ 1479 ierr = PetscMalloc((mbs+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 1480 PetscLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1481 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1482 1483 b->bs = bs; 1484 b->bs2 = bs2; 1485 b->s_nz = 0; 1486 b->s_maxnz = s_nz*bs2; 1487 1488 b->inew = 0; 1489 b->jnew = 0; 1490 b->anew = 0; 1491 b->a2anew = 0; 1492 b->permute = PETSC_FALSE; 1493 PetscFunctionReturn(0); 1494 } 1495 1496 1497 #undef __FUNCT__ 1498 #define __FUNCT__ "MatCreateSeqSBAIJ" 1499 /*@C 1500 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1501 compressed row) format. For good matrix assembly performance the 1502 user should preallocate the matrix storage by setting the parameter nz 1503 (or the array nnz). By setting these parameters accurately, performance 1504 during matrix assembly can be increased by more than a factor of 50. 1505 1506 Collective on MPI_Comm 1507 1508 Input Parameters: 1509 + comm - MPI communicator, set to PETSC_COMM_SELF 1510 . bs - size of block 1511 . m - number of rows, or number of columns 1512 . nz - number of block nonzeros per block row (same for all rows) 1513 - nnz - array containing the number of block nonzeros in the various block rows 1514 (possibly different for each block row) or PETSC_NULL 1515 1516 Output Parameter: 1517 . A - the symmetric matrix 1518 1519 Options Database Keys: 1520 . -mat_no_unroll - uses code that does not unroll the loops in the 1521 block calculations (much slower) 1522 . -mat_block_size - size of the blocks to use 1523 1524 Level: intermediate 1525 1526 Notes: 1527 The block AIJ format is compatible with standard Fortran 77 1528 storage. That is, the stored row and column indices can begin at 1529 either one (as in Fortran) or zero. See the users' manual for details. 1530 1531 Specify the preallocated storage with either nz or nnz (not both). 1532 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1533 allocation. For additional details, see the users manual chapter on 1534 matrices. 1535 1536 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1537 @*/ 1538 int MatCreateSeqSBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1539 { 1540 int ierr; 1541 1542 PetscFunctionBegin; 1543 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1544 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1545 ierr = MatSeqSBAIJSetPreallocation(*A,bs,nz,nnz);CHKERRQ(ierr); 1546 PetscFunctionReturn(0); 1547 } 1548 1549 #undef __FUNCT__ 1550 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1551 int MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1552 { 1553 Mat C; 1554 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1555 int i,len,mbs = a->mbs,nz = a->s_nz,bs2 =a->bs2,ierr; 1556 1557 PetscFunctionBegin; 1558 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1559 1560 *B = 0; 1561 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 1562 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 1563 c = (Mat_SeqSBAIJ*)C->data; 1564 1565 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1566 C->preallocated = PETSC_TRUE; 1567 C->factor = A->factor; 1568 c->row = 0; 1569 c->icol = 0; 1570 c->saved_values = 0; 1571 c->keepzeroedrows = a->keepzeroedrows; 1572 C->assembled = PETSC_TRUE; 1573 1574 c->bs = a->bs; 1575 c->bs2 = a->bs2; 1576 c->mbs = a->mbs; 1577 c->nbs = a->nbs; 1578 1579 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 1580 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 1581 for (i=0; i<mbs; i++) { 1582 c->imax[i] = a->imax[i]; 1583 c->ilen[i] = a->ilen[i]; 1584 } 1585 1586 /* allocate the matrix space */ 1587 c->singlemalloc = PETSC_TRUE; 1588 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1589 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 1590 c->j = (int*)(c->a + nz*bs2); 1591 c->i = c->j + nz; 1592 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));CHKERRQ(ierr); 1593 if (mbs > 0) { 1594 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(int));CHKERRQ(ierr); 1595 if (cpvalues == MAT_COPY_VALUES) { 1596 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1597 } else { 1598 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1599 } 1600 } 1601 1602 PetscLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1603 c->sorted = a->sorted; 1604 c->roworiented = a->roworiented; 1605 c->nonew = a->nonew; 1606 1607 if (a->diag) { 1608 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 1609 PetscLogObjectMemory(C,(mbs+1)*sizeof(int)); 1610 for (i=0; i<mbs; i++) { 1611 c->diag[i] = a->diag[i]; 1612 } 1613 } else c->diag = 0; 1614 c->s_nz = a->s_nz; 1615 c->s_maxnz = a->s_maxnz; 1616 c->solve_work = 0; 1617 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1618 c->mult_work = 0; 1619 *B = C; 1620 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1621 PetscFunctionReturn(0); 1622 } 1623 1624 EXTERN_C_BEGIN 1625 #undef __FUNCT__ 1626 #define __FUNCT__ "MatLoad_SeqSBAIJ" 1627 int MatLoad_SeqSBAIJ(PetscViewer viewer,MatType type,Mat *A) 1628 { 1629 Mat_SeqSBAIJ *a; 1630 Mat B; 1631 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1; 1632 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,*s_browlengths,maskcount; 1633 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1634 int *masked,nmask,tmp,bs2,ishift; 1635 PetscScalar *aa; 1636 MPI_Comm comm = ((PetscObject)viewer)->comm; 1637 1638 PetscFunctionBegin; 1639 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1640 bs2 = bs*bs; 1641 1642 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1643 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1644 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1645 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1646 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1647 M = header[1]; N = header[2]; nz = header[3]; 1648 1649 if (header[3] < 0) { 1650 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 1651 } 1652 1653 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1654 1655 /* 1656 This code adds extra rows to make sure the number of rows is 1657 divisible by the blocksize 1658 */ 1659 mbs = M/bs; 1660 extra_rows = bs - M + bs*(mbs); 1661 if (extra_rows == bs) extra_rows = 0; 1662 else mbs++; 1663 if (extra_rows) { 1664 PetscLogInfo(0,"MatLoad_SeqSBAIJ:Padding loaded matrix to match blocksize\n"); 1665 } 1666 1667 /* read in row lengths */ 1668 ierr = PetscMalloc((M+extra_rows)*sizeof(int),&rowlengths);CHKERRQ(ierr); 1669 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1670 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1671 1672 /* read in column indices */ 1673 ierr = PetscMalloc((nz+extra_rows)*sizeof(int),&jj);CHKERRQ(ierr); 1674 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1675 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1676 1677 /* loop over row lengths determining block row lengths */ 1678 ierr = PetscMalloc(mbs*sizeof(int),&browlengths);CHKERRQ(ierr); 1679 ierr = PetscMalloc(mbs*sizeof(int),&s_browlengths);CHKERRQ(ierr); 1680 ierr = PetscMemzero(s_browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1681 ierr = PetscMalloc(2*mbs*sizeof(int),&mask);CHKERRQ(ierr); 1682 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1683 masked = mask + mbs; 1684 rowcount = 0; nzcount = 0; 1685 for (i=0; i<mbs; i++) { 1686 nmask = 0; 1687 for (j=0; j<bs; j++) { 1688 kmax = rowlengths[rowcount]; 1689 for (k=0; k<kmax; k++) { 1690 tmp = jj[nzcount++]/bs; /* block col. index */ 1691 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 1692 } 1693 rowcount++; 1694 } 1695 s_browlengths[i] += nmask; 1696 browlengths[i] = 2*s_browlengths[i]; 1697 1698 /* zero out the mask elements we set */ 1699 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1700 } 1701 1702 /* create our matrix */ 1703 ierr = MatCreateSeqSBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1704 B = *A; 1705 a = (Mat_SeqSBAIJ*)B->data; 1706 1707 /* set matrix "i" values */ 1708 a->i[0] = 0; 1709 for (i=1; i<= mbs; i++) { 1710 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 1711 a->ilen[i-1] = s_browlengths[i-1]; 1712 } 1713 a->s_nz = 0; 1714 for (i=0; i<mbs; i++) a->s_nz += s_browlengths[i]; 1715 1716 /* read in nonzero values */ 1717 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 1718 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1719 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1720 1721 /* set "a" and "j" values into matrix */ 1722 nzcount = 0; jcount = 0; 1723 for (i=0; i<mbs; i++) { 1724 nzcountb = nzcount; 1725 nmask = 0; 1726 for (j=0; j<bs; j++) { 1727 kmax = rowlengths[i*bs+j]; 1728 for (k=0; k<kmax; k++) { 1729 tmp = jj[nzcount++]/bs; /* block col. index */ 1730 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 1731 } 1732 } 1733 /* sort the masked values */ 1734 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1735 1736 /* set "j" values into matrix */ 1737 maskcount = 1; 1738 for (j=0; j<nmask; j++) { 1739 a->j[jcount++] = masked[j]; 1740 mask[masked[j]] = maskcount++; 1741 } 1742 1743 /* set "a" values into matrix */ 1744 ishift = bs2*a->i[i]; 1745 for (j=0; j<bs; j++) { 1746 kmax = rowlengths[i*bs+j]; 1747 for (k=0; k<kmax; k++) { 1748 tmp = jj[nzcountb]/bs ; /* block col. index */ 1749 if (tmp >= i){ 1750 block = mask[tmp] - 1; 1751 point = jj[nzcountb] - bs*tmp; 1752 idx = ishift + bs2*block + j + bs*point; 1753 a->a[idx] = aa[nzcountb]; 1754 } 1755 nzcountb++; 1756 } 1757 } 1758 /* zero out the mask elements we set */ 1759 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1760 } 1761 if (jcount != a->s_nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 1762 1763 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1764 ierr = PetscFree(browlengths);CHKERRQ(ierr); 1765 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 1766 ierr = PetscFree(aa);CHKERRQ(ierr); 1767 ierr = PetscFree(jj);CHKERRQ(ierr); 1768 ierr = PetscFree(mask);CHKERRQ(ierr); 1769 1770 B->assembled = PETSC_TRUE; 1771 ierr = MatView_Private(B);CHKERRQ(ierr); 1772 PetscFunctionReturn(0); 1773 } 1774 EXTERN_C_END 1775 1776 #undef __FUNCT__ 1777 #define __FUNCT__ "MatRelax_SeqSBAIJ" 1778 int MatRelax_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,int lits,Vec xx) 1779 { 1780 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1781 MatScalar *aa=a->a,*v,*v1; 1782 PetscScalar *x,*b,*t,sum,d; 1783 int m=a->mbs,bs=a->bs,*ai=a->i,*aj=a->j,ierr; 1784 int nz,nz1,*vj,*vj1,i; 1785 1786 PetscFunctionBegin; 1787 its = its*lits; 1788 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %d and local its %d both positive",its,lits); 1789 1790 if (bs > 1) 1791 SETERRQ(PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented"); 1792 1793 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1794 if (xx != bb) { 1795 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 1796 } else { 1797 b = x; 1798 } 1799 1800 ierr = PetscMalloc(m*sizeof(PetscScalar),&t);CHKERRQ(ierr); 1801 1802 if (flag & SOR_ZERO_INITIAL_GUESS) { 1803 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1804 for (i=0; i<m; i++) 1805 t[i] = b[i]; 1806 1807 for (i=0; i<m; i++){ 1808 d = *(aa + ai[i]); /* diag[i] */ 1809 v = aa + ai[i] + 1; 1810 vj = aj + ai[i] + 1; 1811 nz = ai[i+1] - ai[i] - 1; 1812 x[i] = omega*t[i]/d; 1813 while (nz--) t[*vj++] -= x[i]*(*v++); /* update rhs */ 1814 PetscLogFlops(2*nz-1); 1815 } 1816 } 1817 1818 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1819 for (i=0; i<m; i++) 1820 t[i] = b[i]; 1821 1822 for (i=0; i<m-1; i++){ /* update rhs */ 1823 v = aa + ai[i] + 1; 1824 vj = aj + ai[i] + 1; 1825 nz = ai[i+1] - ai[i] - 1; 1826 while (nz--) t[*vj++] -= x[i]*(*v++); 1827 PetscLogFlops(2*nz-1); 1828 } 1829 for (i=m-1; i>=0; i--){ 1830 d = *(aa + ai[i]); 1831 v = aa + ai[i] + 1; 1832 vj = aj + ai[i] + 1; 1833 nz = ai[i+1] - ai[i] - 1; 1834 sum = t[i]; 1835 while (nz--) sum -= x[*vj++]*(*v++); 1836 PetscLogFlops(2*nz-1); 1837 x[i] = (1-omega)*x[i] + omega*sum/d; 1838 } 1839 } 1840 its--; 1841 } 1842 1843 while (its--) { 1844 /* 1845 forward sweep: 1846 for i=0,...,m-1: 1847 sum[i] = (b[i] - U(i,:)x )/d[i]; 1848 x[i] = (1-omega)x[i] + omega*sum[i]; 1849 b = b - x[i]*U^T(i,:); 1850 1851 */ 1852 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1853 for (i=0; i<m; i++) 1854 t[i] = b[i]; 1855 1856 for (i=0; i<m; i++){ 1857 d = *(aa + ai[i]); /* diag[i] */ 1858 v = aa + ai[i] + 1; v1=v; 1859 vj = aj + ai[i] + 1; vj1=vj; 1860 nz = ai[i+1] - ai[i] - 1; nz1=nz; 1861 sum = t[i]; 1862 while (nz1--) sum -= (*v1++)*x[*vj1++]; 1863 x[i] = (1-omega)*x[i] + omega*sum/d; 1864 while (nz--) t[*vj++] -= x[i]*(*v++); 1865 PetscLogFlops(4*nz-2); 1866 } 1867 } 1868 1869 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1870 /* 1871 backward sweep: 1872 b = b - x[i]*U^T(i,:), i=0,...,n-2 1873 for i=m-1,...,0: 1874 sum[i] = (b[i] - U(i,:)x )/d[i]; 1875 x[i] = (1-omega)x[i] + omega*sum[i]; 1876 */ 1877 for (i=0; i<m; i++) 1878 t[i] = b[i]; 1879 1880 for (i=0; i<m-1; i++){ /* update rhs */ 1881 v = aa + ai[i] + 1; 1882 vj = aj + ai[i] + 1; 1883 nz = ai[i+1] - ai[i] - 1; 1884 while (nz--) t[*vj++] -= x[i]*(*v++); 1885 PetscLogFlops(2*nz-1); 1886 } 1887 for (i=m-1; i>=0; i--){ 1888 d = *(aa + ai[i]); 1889 v = aa + ai[i] + 1; 1890 vj = aj + ai[i] + 1; 1891 nz = ai[i+1] - ai[i] - 1; 1892 sum = t[i]; 1893 while (nz--) sum -= x[*vj++]*(*v++); 1894 PetscLogFlops(2*nz-1); 1895 x[i] = (1-omega)*x[i] + omega*sum/d; 1896 } 1897 } 1898 } 1899 1900 ierr = PetscFree(t); CHKERRQ(ierr); 1901 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1902 if (bb != xx) { 1903 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 1904 } 1905 1906 PetscFunctionReturn(0); 1907 } 1908 1909 1910 1911 1912 1913 1914