# This is a BitKeeper generated diff -Nru style patch. # # ChangeSet # 2004/09/13 22:36:02-04:00 victor@tibook500ve.cs.utk.edu # better error handling of lapack; # main routine now returns error # # spai_3.0/lib/spai_error.h # 2004/09/13 22:36:01-04:00 victor@tibook500ve.cs.utk.edu +7 -0 # # spai_3.0/lib/spai_sub.h # 2004/09/13 22:36:01-04:00 victor@tibook500ve.cs.utk.edu +2 -2 # better error handling # # spai_3.0/lib/spai_sub.c # 2004/09/13 22:36:01-04:00 victor@tibook500ve.cs.utk.edu +19 -43 # better error handling # # spai_3.0/lib/spai_error.h # 2004/09/13 22:36:01-04:00 victor@tibook500ve.cs.utk.edu +0 -0 # BitKeeper file /Users/victor/Software/SPAI/orig/spai3-edited/spai_3.0/lib/spai_error.h # # spai_3.0/lib/spai.h # 2004/09/13 22:36:01-04:00 victor@tibook500ve.cs.utk.edu +6 -5 # better error handling # # spai_3.0/lib/spai.c # 2004/09/13 22:36:01-04:00 victor@tibook500ve.cs.utk.edu +34 -19 # better error handling # # spai_3.0/lib/qr.h # 2004/09/13 22:36:01-04:00 victor@tibook500ve.cs.utk.edu +2 -2 # better error handling # # spai_3.0/lib/qr.c # 2004/09/13 22:36:01-04:00 victor@tibook500ve.cs.utk.edu +46 -78 # better error handling # # spai_3.0/lib/mmio.c # 2004/09/13 22:36:01-04:00 victor@tibook500ve.cs.utk.edu +1 -1 # malloc.h removed # # spai_3.0/lib/basics.h # 2004/09/13 22:36:01-04:00 victor@tibook500ve.cs.utk.edu +1 -2 # malloc.h removed # diff -Nru a/spai_3.0/lib/basics.h b/spai_3.0/lib/basics.h --- a/spai_3.0/lib/basics.h 2004-09-13 22:36:38 -04:00 +++ b/spai_3.0/lib/basics.h 2004-09-13 22:36:38 -04:00 @@ -6,11 +6,10 @@ #ifndef __basics_H #define __basics_H +#include #include -#include #include #include -#include #include #ifdef MPI diff -Nru a/spai_3.0/lib/mmio.c b/spai_3.0/lib/mmio.c --- a/spai_3.0/lib/mmio.c 2004-09-13 22:36:38 -04:00 +++ b/spai_3.0/lib/mmio.c 2004-09-13 22:36:38 -04:00 @@ -22,9 +22,9 @@ expressed or implied warranty. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include #include #include -#include #include #include "mmio.h" diff -Nru a/spai_3.0/lib/qr.c b/spai_3.0/lib/qr.c --- a/spai_3.0/lib/qr.c 2004-09-13 22:36:38 -04:00 +++ b/spai_3.0/lib/qr.c 2004-09-13 22:36:38 -04:00 @@ -24,7 +24,7 @@ /**********************************************************************/ -void qr +int qr (matrix *A, int col, int nbq, @@ -54,7 +54,8 @@ if (nbq >= 1) { - multq(Tchar, nbq-1, Ahat,nline_A, block_width); + info = multq(Tchar, nbq-1, Ahat,nline_A, block_width); + if (info) return info; for (j=0, jj=maxapi*dimr, @@ -90,7 +91,7 @@ printf(" n2[nbq]=%d must be <= max=%d\n", n2[nbq],max); printf(".. Try increasing the -mb parameter\n"); - exit(1); + return SPAI_BLAS_ERROR; } TAU_ptr[nbq+1] = TAU_ptr[nbq] + MIN(n2[nbq],nline_A); @@ -101,25 +102,18 @@ } -#ifdef T3D +#if defined(T3D) SGEQRF(&n2[nbq],&nline_A,Z,&max,&TAU[TAU_ptr[nbq]],rw,&max,&info); - if (info) fprintf(stdout,"problem in SGEQRF %d",info); -#endif - -#ifdef ORIGIN - dgeqrf_(&n2[nbq],&nline_A,Z,&max,&TAU[TAU_ptr[nbq]],rw,&max,&info); - if (info) fprintf(stdout,"problem in DGEQRF %d",info); -#endif - -#ifdef INDY +#elif defined(SP2) + dgeqrf(&n2[nbq],&nline_A,Z,&max,&TAU[TAU_ptr[nbq]],rw,&max,&info); +#else dgeqrf_(&n2[nbq],&nline_A,Z,&max,&TAU[TAU_ptr[nbq]],rw,&max,&info); - if (info) fprintf(stdout,"problem in DGEQRF %d",info); #endif -#ifdef SP2 - dgeqrf(&n2[nbq],&nline_A,Z,&max,&TAU[TAU_ptr[nbq]],rw,&max,&info); - if (info) fprintf(stdout,"problem in DGEQRF %d",info); -#endif + if (info) { + fprintf(stdout,"problem in DGEQRF %d",info); + return SPAI_BLAS_ERROR; + } /* C ------- copy Q and R into arrays... */ @@ -182,7 +176,8 @@ else x[jj + i] = 0.0; } - multq(Tchar, nbq, x,nrhs, block_width); + info = multq(Tchar, nbq, x,nrhs, block_width); + if (info) return info; /* compute residual */ fill_zeros(block_width,max,unblocked_len,res); @@ -197,33 +192,28 @@ } } - multq(Nchar, nbq,res,nrhs, block_width); + info = multq(Nchar, nbq,res,nrhs, block_width); + if (info) return info; /* solve upper triangular system, solution in x */ -#ifdef T3D +#if defined(T3D) STRTRS(Uchar_fcd,Nchar_fcd,Nchar_fcd, &dimr,&nrhs,R,&maxapi,x,&max,&info); -#endif - -#ifdef ORIGIN - dtrtrs_(Uchar,Nchar,Nchar, +#elif defined(SP2) + dtrtrs(Uchar,Nchar,Nchar, &dimr,&nrhs,R,&maxapi,x,&max,&info); -#endif - -#ifdef INDY +#else dtrtrs_(Uchar,Nchar,Nchar, &dimr,&nrhs,R,&maxapi,x,&max,&info); #endif -#ifdef SP2 - dtrtrs(Uchar,Nchar,Nchar, - &dimr,&nrhs,R,&maxapi,x,&max,&info); -#endif - - if (info) fprintf(stdout,"problems in DTRTRS %d",info); + if (info) { + fprintf(stdout,"problems in DTRTRS %d",info); + return SPAI_BLAS_ERROR; + } } - + return 0; } /**********************************************************************/ @@ -247,7 +237,7 @@ /**********************************************************************/ -void multq +int multq (char *trans, int nbq, double *A, @@ -279,41 +269,30 @@ /* perform multiplication */ dif = n1[iq+1] - n1[iq]; -#ifdef T3D +#if defined(T3D) SORMQR(Lchar_fcd,Tchar_fcd,&n2[iq],&nline_A,&dif, Qlist[iq], &n2[iq], &TAU[TAU_ptr[iq]], Z,&max,rw,&max,&info); - if (info) fprintf(stdout,"problem in SORMQR\n"); -#endif - -#ifdef ORIGIN - dormqr_(Lchar,Tchar,&n2[iq],&nline_A,&dif, - Qlist[iq], - &n2[iq], - &TAU[TAU_ptr[iq]], +#elif defined(SP2) + dormqr(Lchar,Tchar,&n2[iq],&nline_A,&dif, + Qlist[iq], + &n2[iq], + &TAU[TAU_ptr[iq]], Z,&max,rw,&max,&info); - if (info) fprintf(stdout,"problem in DORMQR\n"); -#endif - -#ifdef INDY +#else dormqr_(Lchar,Tchar,&n2[iq],&nline_A,&dif, Qlist[iq], &n2[iq], &TAU[TAU_ptr[iq]], Z,&max,rw,&max,&info); - if (info) fprintf(stdout,"problem in DORMQR\n"); #endif -#ifdef SP2 - dormqr(Lchar,Tchar,&n2[iq],&nline_A,&dif, - Qlist[iq], - &n2[iq], - &TAU[TAU_ptr[iq]], - Z,&max,rw,&max,&info); - if (info) fprintf(stdout,"problem in DORMQR\n"); -#endif + if (info) { + fprintf(stdout,"problem in DORMQR\n"); + return SPAI_BLAS_ERROR; + } for (j=0, jj=0; jmyid == 0)) printf("\nConstructing scalar SPAI matrix M from A\n"); - M = spai(A, + ierr = spai(A, &M, messages_arg, epsilon_arg, nbsteps_arg, @@ -118,6 +118,8 @@ maxnew_arg, cache_size_arg, verbose_arg); + if (ierr) { + return ierr;} #ifdef MPI MPI_Barrier(A->comm); #endif @@ -138,7 +140,7 @@ if (verbose_arg && (A->myid == 0)) printf("Constructing block SPAI matrix M_B from B\n"); - M_B = spai(B, + ierr = spai(B, &M_B, messages_arg, epsilon_arg, nbsteps_arg, @@ -146,11 +148,15 @@ maxnew_arg, cache_size_arg, verbose_arg); - + if (ierr) { + return ierr;} #ifdef MPI MPI_Barrier(A->comm); #endif + /* VE */ + if (num_bad_cols) printf("SPAI: # bad columns = %d\n",num_bad_cols); + if (verbose_arg) matrix_statistics(M_B,"M_B"); /* Convert M_B to scalar matrix. */ @@ -165,14 +171,15 @@ if (verbose_arg) matrix_statistics(M,"M"); - return(M); + *bspai_matrix = M; + return 0; } /**********************************************************************/ -matrix *spai -(matrix *A, +int spai +(matrix *A, matrix **spai_mat, FILE *messages_arg, /* file for warning messages */ double epsilon_arg, /* tolerance */ int nbsteps_arg, /* max number of "improvement" steps per line */ @@ -183,7 +190,7 @@ int verbose_arg) { matrix *M; - int col; + int col,ierr; int cache_sizes[6]; @@ -255,8 +262,10 @@ */ if (verbose_arg) start_timer(ident_spai); - precompute_column_square_inverses(A); - + ierr = precompute_column_square_inverses(A); + if (ierr) { + return ierr; + } #ifdef MPI MPI_Barrier(A->comm); #endif @@ -271,7 +280,9 @@ } if (col < 0 ) break; - spai_line(A,col,M); + ierr = spai_line(A,col,M); + if (ierr) { + return ierr;} } @@ -301,13 +312,13 @@ if (resplot_fptr) fclose(resplot_fptr); - return(M); - + *spai_mat = M; + return 0; } /**********************************************************************/ -void spai_line +int spai_line (matrix *A, int col, matrix *M) @@ -315,7 +326,7 @@ int s,nbq,nnz,dimr,block_width; double scalar_resnorm,block_resnorm,adjust_epsilon; - int i,index,pe,len; + int i,index,pe,len, ierr; int row_address; int *rptr; @@ -352,7 +363,9 @@ n2[s] = I->slen - dimr; /* compute solution -> x, residual, and update QR */ - qr(A,col,nbq,dimr); + ierr = qr(A,col,nbq,dimr); + if (ierr) { + return ierr;} nnz = J->len; dimr = J->slen; @@ -386,6 +399,7 @@ if (block_resnorm > adjust_epsilon) { num_bad_cols++; + /* if (message) { fprintf(message, "could not meet tol, col=%d resnorm = %le, adjust_epsilon = %le\n", @@ -394,6 +408,7 @@ adjust_epsilon); fflush(message); } + */ } if (resplot_fptr) { @@ -428,7 +443,7 @@ } else break; } - + return 0; } /**********************************************************************/ diff -Nru a/spai_3.0/lib/spai.h b/spai_3.0/lib/spai.h --- a/spai_3.0/lib/spai.h 2004-09-13 22:36:38 -04:00 +++ b/spai_3.0/lib/spai.h 2004-09-13 22:36:38 -04:00 @@ -13,6 +13,7 @@ #include "matrix.h" #include "bicgstab.h" #include "spai_sub.h" +#include "spai_error.h" #include "qr.h" #include "hash_table.h" #include "com_server.h" @@ -89,8 +90,8 @@ extern double **A_addresses; #endif -matrix *bspai -(matrix *, +int bspai +(matrix *, matrix **, FILE *, double, int, @@ -100,8 +101,8 @@ int, int); -matrix *spai -(matrix *, +int spai +(matrix *, matrix **, FILE *, double, int, @@ -110,7 +111,7 @@ int, int); -void spai_line +int spai_line (matrix *, int, matrix *); diff -Nru a/spai_3.0/lib/spai_error.h b/spai_3.0/lib/spai_error.h --- /dev/null Wed Dec 31 16:00:00 196900 +++ b/spai_3.0/lib/spai_error.h 2004-09-13 22:36:38 -04:00 @@ -0,0 +1,7 @@ +#ifndef SPAI_ERROR_H +#define SPAI_ERROR_H + +/* error codes */ +#define SPAI_BLAS_ERROR 1 + +#endif diff -Nru a/spai_3.0/lib/spai_sub.c b/spai_3.0/lib/spai_sub.c --- a/spai_3.0/lib/spai_sub.c 2004-09-13 22:36:38 -04:00 +++ b/spai_3.0/lib/spai_sub.c 2004-09-13 22:36:38 -04:00 @@ -29,7 +29,7 @@ /**********************************************************************/ -void precompute_column_square_inverses +int precompute_column_square_inverses (matrix *A) { int i,j,b,b2,len; @@ -38,7 +38,7 @@ double *work; int local_buf_size,buf_size,index,pe; double *local_si_buffer; - int *buf_sizes,*buf_offsets; + int *buf_sizes,*buf_offsets,ierr; local_buf_size = 0; for (j=0; jmnl; j++) { @@ -56,7 +56,7 @@ n = A->block_sizes[A->my_start_index+j]; Aj_sq_inv = &local_si_buffer[index]; - precompute_column_square_inverse1 + ierr = precompute_column_square_inverse1 (A->lines->A[j], A->lines->ptrs[j], len, @@ -65,6 +65,9 @@ Aj_sq_inv, ipiv, work); + if (ierr) { + return ierr; + } index += (n*n); @@ -117,12 +120,12 @@ free(local_si_buffer); free(ipiv); free(work); - + return 0; } /**********************************************************************/ -void precompute_column_square_inverse1 +int precompute_column_square_inverse1 (double *Aj, int *ptrs, int len, @@ -138,60 +141,33 @@ A,n, Aj_sq_inv); -#ifdef T3D +#if defined(T3D) SGETRF(&n, &n, Aj_sq_inv, &n, ipiv, &info); -#endif - -#ifdef ORIGIN - dgetrf_(&n, &n, Aj_sq_inv, &n, ipiv, &info); -#endif - -#ifdef INDY - dgetrf_(&n, &n, Aj_sq_inv, &n, ipiv, &info); -#endif - -#ifdef SP2 +#elif defined(SP2) dgetrf(&n, &n, Aj_sq_inv, &n, ipiv, &info); +#else + dgetrf_(&n, &n, Aj_sq_inv, &n, ipiv, &info); #endif if (info) { fprintf(stdout,"problem in DGETRF %d\n",info); - exit(1); + return SPAI_BLAS_ERROR; } -#ifdef T3D +#if defined(T3D) DGETRI(&n, Aj_sq_inv, &n, ipiv, work, &n, &info); -#endif - if (info) { - fprintf(stdout,"problem in DGETRI %d\n",info); - exit(1); - } - -#ifdef ORIGIN - dgetri_(&n, Aj_sq_inv, &n, ipiv, work, &n, &info); -#endif - if (info) { - fprintf(stdout,"problem in DGETRI %d\n",info); - exit(1); - } - -#ifdef INDY +#elif defined(SP2) + dgetri(&n, Aj_sq_inv, &n, ipiv, work, &n, &info); +#else dgetri_(&n, Aj_sq_inv, &n, ipiv, work, &n, &info); #endif - if (info) { - fprintf(stdout,"problem in DGETRI %d\n",info); - exit(1); - } -#ifdef SP2 - dgetri(&n, Aj_sq_inv, &n, ipiv, work, &n, &info); -#endif if (info) { fprintf(stdout,"problem in DGETRI %d\n",info); fflush(stdout); - exit(1); + return SPAI_BLAS_ERROR; } - + return 0; } /**********************************************************************/ diff -Nru a/spai_3.0/lib/spai_sub.h b/spai_3.0/lib/spai_sub.h --- a/spai_3.0/lib/spai_sub.h 2004-09-13 22:36:38 -04:00 +++ b/spai_3.0/lib/spai_sub.h 2004-09-13 22:36:38 -04:00 @@ -18,10 +18,10 @@ (index_set *, matrix *); -void precompute_column_square_inverses +int precompute_column_square_inverses (matrix *); -void precompute_column_square_inverse1 +int precompute_column_square_inverse1 (double *, int *, int,