blob: 914ec59ee467a5d1c98cb30294f68eaa5b793b6c [file] [log] [blame]
 /*************************************************************** ****************************************************************** **** Test routine for SGEFA.C **** ****************************************************************** ****************************************************************/ #include #include #include "ge.h" #define BIG 1.0e38 /* Largest representable float. */ #define SMALL 1.0e-45 /* Smallest (in absolute value) representable float. */ #define MATPRINT 7 /* Matricies of size <= MATPRINT are printed out. */ #define VECPRINT 7 /* Vectors of size <= VECPRINT are printed out. */ #ifndef SCALE #define SCALE 1 /* System sizes are scaled by this amount. */ #endif int matgen(); int matvec(); int get_space(); void exit(); main() { register int i, j, k; struct FULL a; /* Storage for rhs, rhs of trans(A)x and solution. */ float *b, *bt, *x, anorm, *col, t; double err, snrm2(); int *ipvt, retval, test_case = 0; void matdump(), ivecdump(), fvecdump(); /* Loop over the test cases. */ /* Initilize matrix. */ while( !matgen( &a, &x, &b, &bt, &ipvt, ++test_case, SCALE ) ) { /* Check that system size is reasonable. */ if( a.rd > MAXCOL || a.cd > MAXCOL ) { fprintf(stderr,"Matrix row dim (%d) or column dim (%d) too big.\n",a.rd,a.cd); exit( 1 ); } /* Calculate the 1 norm of A. */ for( j=0, anorm=0.0; j t ? anorm : t ); } printf("One-Norm(A) ---------- %e.\n", anorm ); /* Test SGEFA. */ retval = sgefa( &a, ipvt ); /* For a successful return from SGEFA test SGESL. */ if( retval ) ;//printf("Zero Column %d found \n", retval ); else { /* Solve system. */ (void)sgesl( &a, ipvt, b, 0 ); if( a.rd <= MATPRINT ) (void) matdump( &a, "FACTORED MATRIX FOLLOWS" ); if( a.rd <= VECPRINT ) { (void)fvecdump( x, a.rd, "True Solution"); (void)fvecdump( b, a.rd, "Solution"); } (void)vexopy( a.rd, b, x, b, 2 ); err = snrm2( a.rd, b, 1 ); //printf(" For Ax=b. Absolute error = %e. Relative error = %e.\n", // err, err/snrm2( a.rd, x, 1 ) ); /* Solve transpose system. */ (void)sgesl( &a, ipvt, bt, 1 ); if( a.rd <= VECPRINT ) { (void)fvecdump( bt, a.rd, "Solution to transposed system"); } (void)vexopy( a.rd, bt, x, bt, 2 ); err = snrm2( a.rd, bt, 1 ); //printf(" For A^Tx=b. Absolute error = %e. Relative error = %e.\n", // err, err/snrm2( a.rd, x, 1 ) ); } } /* End of while loop over test cases. */ return 0; } /* End of MAIN */ int matgen( a, x, b, bt, ipvt, test_case, scale ) struct FULL *a; float **x, **b, **bt; int **ipvt, test_case, scale; /* This routine generates test matrices for the SGE routines. INPUT test_case Switch to type of matrix to generate. 1, 2, 3 => Hilbert slices of various sizes. Note due to the memory allocalion local to this routine test_case=1 must be run before any of the others. 4, 5 => monoelemental test. scale Sizes of systems are scaled according to scale. OUTPUT a The matrix stored in the FULL structure. x Generated solution. b Generated right hand side. bt Generated rhs of trans(A)x. ipvt A pointer to an array of ints. RETURN VALUE 0 => everything went O.K. */ { register int i, j; int n; float *col, tl, tu; char *malloc(); void free(), matdump(), ivecdump(), fvecdump(); if( test_case>1 ) { /* Free up memory used in the last test. */ printf("\n\n**********************************************************************\n"); for( j=0; jrd; j++ ) free( a->pd[j] ); free( *x ); free( *b ); free( *bt ); } /* Determine the test to generate. */ switch( test_case ) { case 1: /* Hilbert slice of various sizes. */ case 2: case 3: n = 3*test_case*scale; /* Set system size. */ a->cd = n; /* Set column and row dimensions. */ a->rd = n; if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */ return( 1 ); printf("Hilbert Slice. Test case %d of size %d.\n", test_case, n ); for( j=0; jpd[j]; i=(j-3) && i<=(j+2) ) *col = 1.0/(float)(i+j+1); } } break; case 4: /* Monoemenental test (NOT SCALED). */ case 5: n=1; a->cd = n; /* Set column and row dimensions. */ a->rd = n; if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */ return( 1 ); printf("Monoelemental. Test case %d of size %d.\n", test_case, n ); *a->pd[0] = ( test_case == 4 ? 3.0 : 0.0 ); break; case 6: /* Tridiagional of various types. */ case 7: case 8: n = 15*scale; a->cd = n; /* Set column and row dimensions. */ a->rd = n; if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */ return( 1 ); printf("Tridiagional. Test case %d of size %d.\n", test_case, n ); tu = 1.0; tl = 1.0; if( test_case == 7 ) tl = 100.0; if( test_case == 8 ) tu = 100.0; for( j=0; jpd[j]; icd = n; /* Set column and row dimensions. */ a->rd = n; if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */ return( 1 ); printf("Rank One. Test case %d of size %d.\n", test_case, n ); for( j=0; jpd[j]; icd = n; /* Set column and row dimensions. */ a->rd = n; if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */ return( 1 ); printf("Zero Column. Test case %d of size %d.\n", test_case, n ); for( j=0; jpd[j]; icd = n; /* Set column and row dimensions. */ a->rd = n; if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */ return( 1 ); printf("Upper Triangular. Test case %d of size %d.\n", test_case, n ); for( j=0; jpd[j]; ij ? 0.0 : (float)(i-j+1) ); break; case 12: /* Lower Triangular. */ n = 6*scale; a->cd = n; /* Set column and row dimensions. */ a->rd = n; if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */ return( 1 ); printf("Lower Triangular. Test case %d of size %d.\n", test_case, n ); for( j=0; jpd[j]; icd = n; /* Set column and row dimensions. */ a->rd = n; if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */ return( 1 ); printf("Near Overflow. Test case %d of size %d. BIG = %e\n", test_case, n, BIG ); tl = (float)(n*n); for( j=0; jpd[j]; ij ? i+1 : j+1 ); /* A number <= 1.0 */ *col = BIG*tu/tl; } break; case 14: /* Near Underflow. */ n = 5*scale; a->cd = n; /* Set column and row dimensions. */ a->rd = n; if( get_space( a, x, b, bt, ipvt ) )/* Get the space needed for vectors. */ return( 1 ); /* Assumes that BIG < 1/SMALL */ printf("Near Underflow. Test case %d of size %d. SMALL = %e\n", test_case, n, 1.0/BIG ); tl = (float)(n*n); for( j=0; jpd[j]; ij ? i+1 : j+1 )/(float)(j+1); /* A number >= 1.0 */ *col = tu*tl/BIG; } break; default: printf("MATGEN: All tests complete.\n"); return( 1 ); break; } /* Generate solution. */ **x = 1.0; if( n>1 ) **bt = 0.0; if( n>2 ) { for( i=2, col=(*x)+2; ird; j++ ) { a->pd[j] = (float *)malloc( a->cd*sizeof( float ) ); if( a->pd[j] == NULL ) { printf("GET_SPACE: Can't get enouph space for matricies...\n"); return( 1 ); } } *x = (float *)malloc( a->cd*sizeof( float ) ); *b = (float *)malloc( a->cd*sizeof( float ) ); *bt = (float *)malloc( a->cd*sizeof( float ) ); *ipvt = (int *)malloc( a->cd*sizeof( int ) ); if( *x == NULL || *b == NULL || *bt == NULL || *ipvt == NULL) { printf("GET_SPACE: Can't get enouph space for vectors...\n"); return( 1 ); } return( 0 ); } int matvec( a, x, b, job) struct FULL *a; float *x, *b; int job; /* This routine calculates b = a*x (if job=0 and b=trans(a)x else) via a column oriented approach. It is most efficient if the matrix is stored by columns. a is a matrix (not its transpose, even when job is nonzero) and x, b are vectors of the appropriate size. RETURNS Nonzero if something goes wrong. */ { register int i, j; float *px, *pb, *col, *row; /* Check input. */ if( (a->cd < 1) || (a->rd < 1) ) return( 1 ); /* Job non-zero => do b = trans(A)x. */ if( job ) { /* Loop over (transposed columns) rows. */ for( i=0, pb=b; ird; i++, pb++ ) { for( j=0, row=a->pd[i], px=x, *pb=0.0; jcd; j++, px++, row++ ) *pb += (*row)*(*px); } return( 0 ); } /* Job zero => do b = Ax. */ /* Loop over columns. */ for( i=0, px=x, pb=b, col=a->pd[0]; ird; i++, pb++, col++) *pb = (*col)*(*px); for( j=1; jcd; j++ ) { for( i=0, px=x+j, pb=b, col=a->pd[j]; ird; i++, pb++, col++) *pb += (*col)*(*px); } return( 0 ); } /* End of MATVEC */ void matdump( a, head ) struct FULL *a; char *head; /* This routine prints out a FULL matrix with headding head. */ { register int i, j; int k, jj, ncolmod, ncolrem; ncolmod = (a->cd)/6; ncolrem = a->cd - ncolmod*6; printf("%s\n", head ); for( i=0; ird; i++) { printf("%3d|", i ); j = 0; for( k=0; k