/************************************************************************/ /* */ /* Program : N A G d i a g I I . h */ /* */ /* Creation : 13/9/2000 */ /* */ /* Class : Library */ /* */ /* Programmer : compiled from NAGLib routines */ /* */ /* Purpose : provide diagonalization routines for other C/C++ */ /* programs. */ /* */ /************************************************************************/ /************************************************************************/ /* */ /* BETA version (works FAPP) */ /* Version II: efs for symm. diag. */ /* */ /************************************************************************/ #include #include #include /*------------------------------------------------------------------------*/ /* NAGLib Types */ /*------------------------------------------------------------------------*/ #define NAG_ERROR_BUF_LEN 512 #define NAG_FILE_LEN 80 #define NAGERR_DEFAULT (NagError*)0 #define COMPLEX typedef char Boolean; typedef struct { double re,im; } Complex; typedef struct { /* Nag Error Structure */ int code; /* Out: Error Code */ Boolean print; /* In: print? yes/no */ char message[NAG_ERROR_BUF_LEN]; /* InOut: Error message */ void (*handler)(char*,int,char*);/* In: Error handling fctn */ long errnum; /* Value for some errors */ } NagError; /*------------------------------------------------------------------------*/ /* Subroutines provided by this code (see User Interface section) */ /*------------------------------------------------------------------------*/ void diagonalize_symmetric(double **matrix, int n, double **evs, double **efs); void diagonalize_asymmetric(double **matrix, int n, Complex **evs, Complex **efs); void matrix2MathematicaFile(double *mat, int n, char* name); void matrix2SparseFile(double *mat, int n, char* name); /*------------------------------------------------------------------------*/ /* nrutil Prototypes */ /*------------------------------------------------------------------------*/ double *dvector(long nl, long nh); /* allocate a double vector with subscript range v[nl..nh] */ Complex *Cvector(long nl, long nh); /* allocate a Complex vector with subscript range v[nl..nh] */ long *Intvector(long nl, long nh); /* allocate a long vector with subscript range v[nl..nh] */ void free_dvector(double *v, long nl, long nh); /* free a double vector allocated with dvector() */ void free_Cvector(Complex *v, long nl, long nh); /* free a Complex vector allocated with Cvector() */ /*------------------------------------------------------------------------*/ /* NAGLib Prototypes */ /*------------------------------------------------------------------------*/ void f02abc(long n, double *a, long tda, double *r, double *v, long tdv, NagError *fail); /* diagonalize real symmetric matrix */ void f02agc(long n, double a[], long tda, Complex r[], Complex v[], long tdv, long iter[], NagError *fail); /* diagonalize real asymmetric matrix */ /*------------------------------------------------------------------------*/ /* User Interface */ /*------------------------------------------------------------------------*/ void diagonalize_symmetric(double **matrix, int n, double **evs, double **efs) /*------------------------------------------------------------------------*/ /* diagonalizes quadratic symmetric matrix of dimension . */ /* On exit: eigenvalues in , eigenvectors as columns in */ /* Usage: */ /* ------ double *Mat, *eigenvalues, *eigenfunctions; */ /* Mat = dvector(0,dim*dim - 1); */ /* ... */ /* diagonalize_symmetric(&Mat,dim,&eigenvalues,&eigenfunctions); */ /* */ /*------------------------------------------------------------------------*/ { int IA_nag,N_nag,IV_nag; /* variables for NAGLib */ double *e,*v,*help; NagError IFAIL_nag; IA_nag=IV_nag=N_nag=n; *evs = dvector(0,n-1); e = dvector(0,n-1); *efs = dvector(0,n*n-1); f02abc((long)n, *matrix, (long)n, *evs, *efs, (long)n, &IFAIL_nag); if (IFAIL_nag.code!=0) printf("Error in : code=%1d\n",IFAIL_nag.code); else printf("In II: diagonalization successful.\n"); free_dvector(e,0,n-1); } void diagonalize_asymmetric(double **matrix, int n, Complex **evs, Complex **efs) /*-------------------------------------------------------------------------*/ /* diagonalizes quadratic asymmetric matrix of dimension . */ /* On exit: eigenvalues in , eigenvectors as columns in */ /* Usage: */ /* ------ double *Mat; */ /* Complex *eigenvalues, *eigenfunctions; */ /* Mat = dvector(0,dim*dim - 1); */ /* ... */ /* diagonalize_asymmetric(&Mat,dim,&eigenvalues,&eigenfunctions); */ /* */ /*-------------------------------------------------------------------------*/ { int IA_nag,N_nag,IV_nag; /* variables for NAGLib */ long *intger; /* variable for unsym. NAGLIB diagonalization */ NagError IFAIL_nag; IA_nag=IV_nag=N_nag=n; intger=Intvector(0,n*n-1); *evs =Cvector(0,n-1); *efs =Cvector(0,n*n-1); f02agc((long)n, *matrix, (long)n, *evs, *efs, (long)n, intger, &IFAIL_nag); if (IFAIL_nag.code!=0) printf("Error in : code=%1d\n",IFAIL_nag.code); else printf("In II: diagonalization successful.\n"); } void matrix2MathematicaFile(double *mat, int n, char* name) /*-------------------------------------------------------------------------*/ /* Write matrix of dim. to file in Mathematica notation. */ /* Executing the file in Mathematica will produce eigenvalues of . */ /* Usage: */ /* ------ double *Mat; */ /* ... */ /* matrix2MathematicaFile(Mat,dimension,"Matrix.mat"); */ /* */ /*-------------------------------------------------------------------------*/ { long i,j,index=0; FILE *fp; fp=fopen(name,"w"); fprintf(fp,"Clear[matrix];\n\n"); fprintf(fp,"matrix={"); for (i=0; i of dim. to file in sparse matrix notation.*/ /* Usage: */ /* ------ double *Mat; */ /* ... */ /* matrix2SparseFile(Mat,dimension,"Matrix.mat"); */ /* */ /*-------------------------------------------------------------------------*/ { long i,j,index=0; FILE *fp; fp=fopen(name,"w"); for (i=1; i<=n; ++i) { fprintf(fp,"%6ld 0.0\n",-i); for (j=1; j<=n; ++j) { if (fabs(mat[index])>FLT_EPSILON) fprintf(fp,"%6ld %18.12f\n",j,mat[index]); ++index; } } fclose(fp); } /*------------------------------------------------------------------------*/ /* End User Interface */ /*------------------------------------------------------------------------*/ #if !defined(NumRec) #include "allocation.h" #endif /*------------------------------------------------------------------------*/ /* NAGLib Routines */ /*------------------------------------------------------------------------*/ /*#include "nag.h"*/ /* * * Copyright 1990 Numerical Algorithms Group * * General include file for all NAG C software. * * Malcolm Cohen, NAG Ltd., Oxford, U.K., 1990. * * Mark 1, 1990. * * Mark 2, revised * The defined types and error names have been placed into * separate files, nag_types.h and nag_errlist.h, respectively. * SPD, October 91. * Reworked handling of implementation specifics. AJS feb92. * * Purpose : To define a macro that is implementation specific so that * variant code may be easily inserted by implementors and carried * foreward into future releases. * A number of other macros are also defined to indicate the * compliance (or otherwise) to particular common characteristics. * These macros are always set to their default state in the * introduction. There is a later section where the implementation- * specific macro controls changes to the defaults. */ #ifndef NAG_H #define NAG_H /* The type of system should be defined here. It is defined by the "system" and "compiler" components of the NAG product code for the library (c.f. a00aact.c). Some examples are : APOA_NAG_IMP for Apollo SR 10.3 IBPA_NAG_IMP for the microsoft PC compiler IBPC_NAG_IMP for the Zortech PC compiler IBPE_NAG_IMP for the Borland PC compiler SU3A_NAG_IMP for Sun 3 These macros are used in preprocessor directives in other parts of the code to introduce the implementation specific variations. It is only anticipated that variations will be required in nag.h , nag_stddef.h , nag_stdlib.h , nag_string.h and x05aact.c. */ /*#define APOA_NAG_IMP Apollo SR 10.3 */ /*#define D31A_NAG_IMP DEC RISC Ultrix */ /*#define DVVA_NAG_IMP DEC VAX/VMS */ /*#define DVVG_NAG_IMP DEC VAX/VMS G_float option */ /*#define IBPA_NAG_IMP MS/PC-DOS PCs , Microsoft compiler */ /*#define IBPC_NAG_IMP MS/PC-DOS PCs , Zortech compiler */ /*#define IBPE_NAG_IMP MS/PC-DOS & MS/WINDOWS PCs , Borland compiler */ /*#define SU3A_NAG_IMP Sun 3 bundled C compiler */ /*#define SU4A_NAG_IMP Sun 4 (Sparc) */ /*#define SG4A_NAG_IMP Silicon Graphics IRIX */ #define DAUA_NAG_IMP DEC Alpha AXP OSF/1 /* The following 13 sections define the default settings for macros concerned with common characteristics that are shared by many compilers. These settings are listed below and a macro is defined for each setting. If for a particular compiler, a setting is not appropriate, then it should be corrected later, in the implementation specific section of this file by undefining the corresponding macro. For example, we assume that most C compilers support function prototypes, hence NAG_PROTO is defined. If the compiler happens to be the bundled SUN C compiler which does not support function prototypes, then this macro is undefined, see below. */ /* NAG_PROTO This is used to indicate that the compiler supports function prototypes as specified in the ANSI standard. For the older K&R compilers which do not support function prototypes undefine NAG_PROTO */ #define NAG_PROTO /* NAG_FABS This macro is defined on the assumption that the function fabs is available in the math library. If this is not so and the compiler does not offer anything better then use the macro ABS which is defined later in this file by undefining NAG_FABS. */ #define NAG_FABS /* NAG_EXIT The default action for p01acc is to terminate the program by the use of exit. On some systems it may be more appropriate to terminate with abort, perhaps causing a core dump. However, any parent process (or script) may wish to continue after testing the exit status of the program. If it is approriate to use abort rather than exit then undefine NAG_EXIT. */ #define NAG_EXIT /* NAG_YES_STDLIB Most compilers provide the include file stdlib.h containining declarations for the functions abort , atoi , free , exit , atol , calloc , malloc , realloc , alloca. If this is not the case then undefine stdlib.h */ #define NAG_YES_STDLIB /* NAG_YES_STRING Most compilers provide the include file string.h, however some older BSD derived, non-ANSI compilers provide strings.h instead. If strings.h is provided then undefine NAG_YES_STRING */ #define NAG_YES_STRING /* NAG_YES_STDDEF Most compilers provide the include file stddef.h containining definitions for ptrdiff_t and size_t. If these definitions are not provided then undefine NAG_YES_STDDEF. */ #define NAG_YES_STDDEF /* NAG_VOID_STAR Newer ANSI based compilers provide the generic pointer void *. If the target system does not allow the type void *, then undefine NAG_VOID_STAR. */ #define NAG_VOID_STAR /* NAG_IND_STR Newer ANSI based compilers provide the string functions strchr and strrchr. Some older compilers provide index instead of strchr and rindex instead of strrchr. If this is the case with the target system then undefine NAG_IND_STR. */ #define NAG_IND_STR /* NAG_NO_INC_MEMORY_DOT_H The functions memcmp , memcpy are normally declared string(s).h. If on target system they are declared in memory.h then undefine NAG_NO_INC_MEMORY_DOT_H */ #define NAG_NO_INC_MEMORY_DOT_H /* NAG_NO_BC_MEM_FUN The functions memcpy and memcmp are meant to manipulate objects as character arrays. Some BSD derived systems provide bcopy and bcmp as either additional to or instead of memcpy and memcmp. If it is felt that it is desirable to use bcopy and bcmp then undefine NAG_NO_BC_MEM_FUN */ #define NAG_NO_BC_MEM_FUN /* NAG_ARRAY_SIZE In the stringent test program for c06fuc, a number of large arrays are being used. This might prove too large, especially for the PC type machines. In this case undefine NAG_ARRAY_SIZE. */ #define NAG_ARRAY_SIZE /* NAG_TIME_DOT_H x05aac function is used to set up a non-repeatable seed for the random number generator function g05ccc. This is achieved by setting the seed to a value derived from the system time by calling an appropriate function. Ansi derived compilers call the function time with the argument of type time_t and include files time.h and sys/types.h. If this is not the case undefine NAG_TIME_DOT_H. System specific code may then have to be inserted in x05aact.c */ #define NAG_TIME_DOT_H /* NAG_YES_MALLOC The C dynamic memory allocator, malloc returns size_t whose maximum value may not be large enough for practical applications, e.g the Microsoft C compiler whose size_t has the maximum value of 65536. In this case the compiler may provide an alternative function, such as halloc. See nag_stdlib.h for the definition of halloc. */ #define NAG_YES_MALLOC /* NAG_LINT The Apollo version of lint is not ANSI compatible. For lint to work successfully, it seems necessary to assume a 'pure' bsd system. */ /* NAG_CHAR_BOOL Most compilers seem to be happy to accept typedef char Boolean; If yours insists on defining typedef signed char Boolean; then undef NAG_CHAR_BOOL */ #define NAG_CHAR_BOOL /* The following sections override the default settings of the above 14 macros, where it is necessary for a particular implementation. */ #ifdef IBPA_NAG_IMP #undef NAG_EXIT #undef NAG_ARRAY_SIZE #undef NAG_YES_MALLOC #else #ifdef IBPE_NAG_IMP #undef NAG_ARRAY_SIZE #undef NAG_YES_MALLOC extern unsigned _floatconvert; #pragma extref _floatconvert #else #ifdef SU3A_NAG_IMP #undef NAG_PROTO #undef NAG_EXIT #undef NAG_NO_BC_MEM_FUN #undef NAG_IND_STR #else #ifdef SU4A_NAG_IMP #undef NAG_PROTO #undef NAG_EXIT #undef NAG_NO_BC_MEM_FUN #undef NAG_IND_STR #else #ifdef DVVA_NAG_IMP #undef NAG_TIME_DOT_H #else #ifdef NAG_LINT #undef NAG_PROTO #undef NAG_VOID_STAR #undef NAG_YES_STRING #undef NAG_YES_STDLIB #undef __STDC__ #else #ifdef SG4A_NAG_IMP #undef NAG_CHAR_BOOL #endif #endif #endif #endif #endif #endif #endif /* ****** End of anticipated implementation specifics ***** */ #ifndef NULLFN #define NULLFN (void(*)())0 #endif #ifndef NULLDFN #define NULLDFN (double(*)())0 #endif #ifndef TRUE #define TRUE 1 #define FALSE 0 #endif #define NAGERR_DEFAULT (NagError*)0 #define NAGUSER_DEFAULT (Nag_User *)0 #define NAGCOMM_NULL (Nag_Comm *)0 #define NAGMESG_DEFAULT (Nag_Mesg *)0 #define E01_DEFAULT (Nag_E01_Opt *)0 #define E04_DEFAULT (Nag_E04_Opt *)0 #define G13_DEFAULT (Nag_G13_Opt *)0 #define TRANSF_ORDER_DEFAULT (Nag_TransfOrder *)0 /* INIT_FAIL initialises a fail structure with print set to FALSE */ #define INIT_FAIL(E) (E).code = NE_NOERROR, (E).print = 0, \ (E).handler = 0, (E).errnum = 0 /* SET_FAIL initialises a fail structure with print set to TRUE */ #define SET_FAIL(E) (E).code = NE_NOERROR, (E).print = TRUE, \ (E).handler = 0, (E).errnum = 0 #define INIT_MESG(M) (M).code = NM_NO_MESG, (M).print = FALSE, \ (M).res_file = (char *)0, (M).name = (char *)0, (M).opt_name = (char *)0, \ (M).print_mesg = 0, (M).init_mesg1 = IDUMMY, (M).init_mesg2 = INIT2DUMMY #define INIT_STREAM(S) (S).set = FALSE, (S).st_out = FALSE, \ (S).st_err = FALSE, (S).st_in = FALSE, (S).chapter = (char *)0, \ (S).error = NE_NOERROR /* d02rac use, define two or four functions as null pointers in the d02rac * call. Saves us having to design a limited option setting facility for * this routine */ #define NULL_4_FUN NULLFN, NULLFN, NULLFN, NULLFN #define NULL_2_FUN NULLFN, NULLFN /* Magic numbers to denote initialisation state */ #define RDUMMY -11111.0 #define IDUMMY -11111 #define INIT2DUMMY -23456 #define Vprintf (void)printf #define Vfprintf (void)fprintf #define Vsprintf (void)sprintf #define Vscanf (void)scanf #define Vfscanf (void)fscanf #define Vstrcpy (void)strcpy #define ABS(x) (x>=0 ? x : -(x)) /* Absolute value */ #ifdef NAG_FABS #define FABS(x) fabs(x) #else #define FABS(x) ABS(x) #endif #define SIGN(x,y) (y>0 ? ABS(x) : -ABS(x)) /* Sign transfer */ #define MAX(x,y) (x>=y ? x : y) /* Maximum of two arguments */ #define MIN(x,y) (x * * Copyright 1990 Numerical Algorithms Group Ltd. * * Include file for NAG C Library * * Mark 1, 1990 * * Mark 2, revised. * Revised handling of implementation specifics. AJS jan92. * * Purpose : Ensures that typedefs for size_t and ptrdiff_t are included. */ #ifdef NAG_YES_STDDEF /*#include "stddef.h"*/ #else #ifndef SIZE_T #define SIZE_T typedef unsigned size_t; #endif /* not SIZE_T */ #ifndef PTRDIFF_T #define PTRDIFF_T typedef int ptrdiff_t; /* result of subtracting two pointers */ #endif /* not PTRDIFF_T */ #endif #endif /* not NAG_STDDEF */ /*------------------------------------------------------------------*/ /* New Header */ /*------------------------------------------------------------------*/ #ifndef NAG_STDLIB #define NAG_STDLIB /* * * Copyright 1989 Numerical Algorithms Group * * NAG routine interface to ANSI stdlib routines. * * Malcolm Cohen, 1989, NAG Ltd, Oxford, U.K. * * Mark 1, 1990. * Mark 1a revised. IER-864 (Oct 1990). */ #ifdef NAG_YES_STDLIB #include #ifndef NAG_STDDEF #define NAG_STDDEF /* * * Copyright 1990 Numerical Algorithms Group Ltd. * * Include file for NAG C Library * * Mark 1, 1990 * * Mark 2, revised. * Revised handling of implementation specifics. AJS jan92. * * Purpose : Ensures that typedefs for size_t and ptrdiff_t are included. */ #ifdef NAG_YES_STDDEF /*#include "stddef.h"*/ #else #ifndef SIZE_T #define SIZE_T typedef unsigned size_t; #endif /* not SIZE_T */ #ifndef PTRDIFF_T #define PTRDIFF_T typedef int ptrdiff_t; /* result of subtracting two pointers */ #endif /* not PTRDIFF_T */ #endif #endif /* not NAG_STDDEF */ #else #include "nag_stddef.h" #ifdef NAG_PROTO extern void abort(void); extern void exit(int status); extern int atoi(const char *s); extern long atol(const char *s); extern void *calloc(size_t nobj, size_t size); extern void *malloc(size_t size); extern void far *farmalloc(unsigned long size); extern void *realloc(void *p, size_t size); extern void far *farrealloc(void far *p, unsigned long size); extern void free(void *p); extern void farfree(void far *p); extern char *alloca(int size); #else extern int abort(),atoi(),free(),farfree(),exit(); extern long atol(); extern char *calloc(),*malloc(),*realloc(),*alloca(); /*extern char far *farmalloc(); extern char far *farrealloc(); */ #endif #endif #ifdef NAG_YES_MALLOC #define NAG_ALLOC(n,type) (type *)malloc((size_t)(n)*sizeof(type)) #define NAG_REALLOC(pointer,n,type) (type *)realloc((Pointer)(pointer), \ (size_t)(n)*sizeof(type)) #define NAG_FREE(x) free((Pointer)x) #else #ifdef IBPA_NAG_IMP #include #define NAG_ALLOC(n,type) (type *)halloc((long)(n), sizeof(type)) #define NAG_REALLOC(pointer,n,type) (type *)realloc((Pointer)(pointer), \ (size_t)(n)*sizeof(type)) #define NAG_FREE(x) hfree((Pointer)x) #else #ifdef IBPE_NAG_IMP #include #include #define NAG_CALLOC(n,type) (type huge *)farcalloc((long)(n), sizeof(type)) #define NAG_ALLOC(n,type) (type far *)farmalloc((long)(n)*sizeof(type)) #define NAG_HALLOC(n,type) (type huge *)farmalloc((long)(n)*sizeof(type)) #define NAG_REALLOC(pointer,n,type) (type far *)farrealloc((Pointer)(pointer),\ (long)(n)*sizeof(type)) #define NAG_FREE(x) farfree((Pointer)x) #endif #endif #endif #ifndef EXIT_SUCCESS #define EXIT_SUCCESS 0 #endif #ifndef EXIT_FAILURE #define EXIT_FAILURE 1 #endif #endif /* not NAG_STDLIB */ /*------------------------------------------------------------------*/ /* New Header */ /*------------------------------------------------------------------*/ /*#include "nagf01.h"*/ /*------------------------------------------------------------------*/ /* New Header */ /*------------------------------------------------------------------*/ /*#include "nagf02.h"*/ /*------------------------------------------------------------------*/ /* New Header */ /*------------------------------------------------------------------*/ /*#include "nagf06.h"*/ #ifdef NAG_PROTO extern void f06aaz(char *srname, long info); #else extern void f06aaz(); #endif /*------------------------------------------------------------------*/ /* New Header */ /*------------------------------------------------------------------*/ #ifndef NAGX02 #define NAGX02 /* * * Copyright 1991 Numerical Algorithms Group. * * Header file for the x02 chapter of the NAG C Library. * Definitions of machine-specific constants. * * This header file was generated by getcons. * * Mark 2, 1991 * */ #define X02AHC x02ahc() #define X02AJC x02ajc() #define X02AKC x02akc() #define X02ALC x02alc() #define X02AMC x02amc() #ifdef NAG_PROTO extern double x02ahc(void); extern double x02ajc(void); extern double x02akc(void); extern double x02alc(void); extern double x02amc(void); #else extern double x02ahc(); extern double x02ajc(); extern double x02akc(); extern double x02alc(); extern double x02amc(); #endif #define X02BBC 9223372036854775807 #define X02BEC 15 #define X02BHC 2 #define X02BJC 53 #define X02BKC -1021 #define X02BLC 1024 #define X02DAC FALSE #define X02DJC TRUE #endif /* not NAGX02 */ /*------------------------------------------------------------------*/ /* New Header */ /*------------------------------------------------------------------*/ #ifndef NAGP01 #define NAGP01 /* * * Copyright 1996 Numerical Algorithms Group * * Include file for NAG C Library p01 Chapter * * Mark 4 re-issue, 1996. */ #ifdef NAG_PROTO extern void p01acc(char *string, int code, char *name, NagError *fail); #else extern void p01acc(); #endif #endif /* not NAGP01 */ /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f06jlct.c * * Copyright 1990 Numerical Algorithms Group. * * NAG C Library. * * Derived from: NAG Fortran Library F06JLF(Last revision Mark 12) * * C version: Ian Hounam, August 1989, NAG Ltd Oxford, U.K. * * Mark 1, 1990 * Mark 2 Global Revision 1991 * Nag/ removed from include file names. SPD. */ #ifdef NAG_PROTO long f06jlc(long n, double *x, long incx) #else long f06jlc(n, x, incx) long n, incx; double *x; #endif { double xmax, tmp; long i, ix, imax; if (n > 0) { imax = 1; if (n > 1) { xmax = ABS(x[0]); ix = 1; for (i = 2; i <= n; ++i) { ix += incx; if (xmax < (tmp = ABS(x[ix-1]))) { xmax = tmp; imax = i; } } } } else { imax = 0; } return imax; } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f06edct.c * * Copyright 1990 Numerical Algorithms Group. * * NAG C Library. * * Derived from: NAG Fortran Library F06EDF(Last revision Mark 12) * * C version: Ian Hounam, August 1989, NAG Ltd Oxford, U.K. * * Mark 1, 1990 * Mark 2 Global Revision 1991 * Nag/ removed from include file names. SPD. */ #ifdef NAG_PROTO void f06edc(long n, double alpha, double *x, long incx) #else void f06edc(n, alpha, x, incx) long n, incx; double alpha, *x; #endif { long ix, mx; if (n > 0) { if (alpha == 0.0) { mx = (n - 1) * incx; for (ix = 0; incx < 0 ? ix >= mx : ix <= mx; ix += incx) x[ix] = 0.0; } else if (alpha == -1.0) { mx = (n - 1) * incx; for (ix = 0; incx < 0 ? ix >= mx : ix <= mx; ix += incx) x[ix] = -x[ix]; } else if (alpha != 1.0) { mx = (n - 1) * incx; for (ix = 0; incx < 0 ? ix >= mx : ix <= mx; ix += incx) { x[ix] = alpha * x[ix]; } } } } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f06eact.c * * Copyright 1990 Numerical Algorithms Group. * * NAG C Library. * * Derived from: NAG Fortran Library F06EAF(Last revision Mark 12) * * C version: Ian Hounam, August 1989, NAG Ltd Oxford, U.K. * * Mark 1, 1990 * Mark 2 Global Revision 1991 * Nag/ removed from include file names. SPD. */ #ifdef NAG_PROTO double f06eac(long n, double *x, long incx, double *y, long incy) #else double f06eac(n, x, incx, y, incy) long n, incx, incy; double *x, *y; #endif { long i; double sum; sum = 0.0; if (n > 0) { if (incy < 0) y += - (n - 1) * incy; if (incx < 0) x += - (n - 1) * incx; if (incx==incy && incx==1) for (i=n; --i>=0; sum += *x++ * *y++); else for (i=n; --i>=0; (x+=incx),(y+=incy)) sum += *x * *y; } return sum; } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f06prct.c * * Copyright 1994 Numerical Algorithms Group. * * NAG C Library. * * dsyr2 level 2 C Blas. * * Author : George Levy, NAG Ltd., Oxford, U.K. * * Mark 3 re-issue, 1994. * */ #ifdef NAG_PROTO void f06prc(MatrixTriangle UpperLower, long n, double alpha, const double x[], long incx, const double y[], long incy, double a[], long tda) #else void f06prc( UpperLower, n, alpha, x, incx, y, incy, a, tda) MatrixTriangle UpperLower; long n; double alpha; double x[]; long incx; double y[]; long incy; double a[]; long tda; #endif { /* Local variables */ long info; double temp1, temp2; long i, j; long ix, iy, jx, jy, kx, ky; #define A(I,J) a[(J)+(I)*tda] /* Purpose * ======= * * dsyr2 performs the symmetric rank 2 operation * * A := ALPHA*X*Y' + ALPHA*Y*X' + A. * * * where ALPHA is a real scalar, X and Y are n element vectors and A is an * n by n symmetric matrix. * * Parameters * ========== * * uplo MatrixTriangle - specifies whether the matrix is an upper or lower * triangular matrix as follows: * uplo = UpperTriangle A is stored in upper triangular matrix. * uplo = LowerTriangle A is stored in lower triangular matrix. * * * n long - specifies the order of the matrix A. * (n must be at least zero). * * * alpha double - specifies the scalar ALPHA. * * * x const double [] - x is a pointer to a vector of dimension at least * (1+(n-1)*ABS(incx)). The first element of the vector is contained in x[0], * and the (i+1)th in x[i]. The incremented array x must contain the n * element vector X. * * * incx long - specifies the increment for the elements of x. * (incx must not be zero). * * * y const double [] - y is a pointer to a vector of dimension at least * (1+(n-1)*ABS(incy)). The first element of the vector is contained in y[0], * and the (i+1)th in y[i]. The incremented array x must contain the n * element vector Y. * * * incy long - specifies the increment for the elements of y. * (incy must not be zero). * * * a double [] - a[0] contains A(0, 0), the element in the first * row and first column of matrix A. The mapping used from a to A is * A(I, J) = a[(J) + (I) * tda ] where I is the row index and J is the column * index. * With uplo = UpperTriangle, the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * With uplo = LowerTriangle, the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * * * tda long - specifies the trailing dimension used in the definition * of matrix A. (tda must be at least MAX(1, n)). * */ info = 0; if (UpperLower!=UpperTriangle && UpperLower!=LowerTriangle) info = 1; else if (n < 0) info = 2; else if (incx == 0) info = 5; else if (incy == 0) info = 7; else if ( tda < MAX(1, n)) info = 9; if (info != 0) { f06aaz("f06prc/dsyr2", info); return; } /* Quick return if possible. */ if (n == 0 || alpha == 0.0) return; /* Set up the start points in x and y if the increments are not both unity. */ if (incx != 1 || incy != 1) { if (incx > 0) kx = 1; else kx = 1 - (n - 1) * incx; if (incy > 0) ky = 1; else ky = 1 - (n - 1) * incy; jx = kx - 1; jy = ky - 1; } /* Start the operations. In this version the elements of A are accessed sequentially with one pass through the triangular part of A. */ if (UpperLower==UpperTriangle) { /* Form A when A is stored in the upper triangle. */ if (incx == 1 && incy == 1) { for (i = 0; i < n; ++i) { if (x[i] != 0.0 || y[i] != 0.0) { temp1 = alpha * x[i]; temp2 = alpha * y[i]; for (j = i; j < n; ++j) A(i, j) += temp1 * y[j] + temp2 * x[j]; } } } else { ix = kx - 1; iy = ky - 1; for (i = 0; i < n; ++i) { if (x[ix] != 0.0 || y[iy] != 0.0) { temp1 = alpha * x[ix]; temp2 = alpha * y[iy]; jx = ix; jy = iy; for (j = i; j < n; ++j) { A(i, j) += temp1 * y[jy] + temp2 * x[jx]; jx += incx; jy += incy; } } ix += incx; iy += incy; } } } else { /* Form A when A is stored in the lower triangle. */ if (incx == 1 && incy == 1) { for (i = 0; i < n; ++i) { if (x[i] != 0.0 || y[i] != 0.0) { temp1 = alpha * x[i]; temp2 = alpha * y[i]; for (j = 0; j <= i; ++j) A(i, j) += temp1 * y[j] + temp2 * x[j]; } } } else { ix = kx - 1; iy = ky - 1; for (i = 0; i < n; ++i) { if (x[ix] != 0.0 || y[iy] != 0.0) { temp1 = alpha * x[ix]; temp2 = alpha * y[iy]; jx = kx - 1; jy = ky - 1; for (j = 0; j <= i; ++j) { A(i, j) += temp1 * y[jy] + temp2 * x[jx]; jx += incx; jy += incy; } } ix += incx; iy += incy; } } } } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f06pmct.c * * Copyright 1994 Numerical Algorithms Group. * * NAG C Library. * * dger level 2 C Blas. * * Author : George Levy, NAG Ltd., Oxford, U.K. * * Mark 3 re-issue, 1994. * */ #ifdef NAG_PROTO void f06pmc(long m, long n, double alpha, const double x[], long incx, const double y[], long incy, double a[], long tda) #else void f06pmc(m, n, alpha, x, incx, y, incy, a, tda) long m, n; double alpha; double x[]; long incx; double y[]; long incy; double a[]; long tda; #endif { /* Local variables */ long info; double temp; long i, j; long ix, jy, kx, ky; #define A(I,J) a[(J)+(I)*tda] /* Purpose * ======= * * dger performs the rank 1 operation * * A := ALPHA*X*Y' + A, * * where ALPHA is a scalar, X is an m element vector, Y is an n * element vector and A is an m by n matrix. * * Parameters * ========== * * m long - specifies the number of rows of the matrix A. * (m must be at least zero). * * * n long - specifies the number of columns of the matrix A. * (n must be at least zero). * * * alpha double - specifies the scalar ALPHA. * * * x const double [] - x is a pointer to a vector of dimension at least * (1+(m-1)*ABS(incx)). The first element of the vector is x[0], and the * (i+1)th is x[i]. The incremented array x must contain the m element * vector X. * * * incx long - specifies the increment for the elements of x. * (incx must not be zero). * * * y const double [] - y is a pointer to a vector of dimension at least * (1+ (n-1)*ABS(incy)). The first element of the vector is y[0], and the * (i+1)th is y[i]. The incremented array y must contain the n * element vector Y. * * * incy long - specifies the increment for the elements of y. * (incy must not be zero). * * * a double [] - a[0] contains A(0, 0), the element in the first row * and first column of matrix A. The mapping used from a to A is * A(I, J) = a[(J) + (I) * tda] where I is the row index, with range * I = 0 .. m - 1, and J is the column index, with range J = 0 .. n - 1. * * * tda long - specifies the trailing dimension used in the definition * of Matrix A. (tda must be at least MAX(1, n)). * */ info = 0; if (m < 0) info = 1; else if (n < 0) info = 2; else if (incx == 0) info = 5; else if (incy == 0) info = 7; else if (tda < MAX(1, n)) info = 9; if (info != 0) { f06aaz("f06pmc/dger", info); return; } /* Quick return if possible. */ if (m == 0 || n == 0 || alpha == 0.0) return; /* Start the operations. In this version the elements of A are accessed sequentially with one pass through A. */ if (incx > 0) kx = 0; else kx = 1 - (m - 1) * incx - 1; if (incy == 1) { ix = kx; for (i = 0; i < m; ++i) { if (x[ix] != 0.0) { temp = alpha * x[ix]; for (j = 0; j < n; ++j) { A(i, j) += temp * y[j]; } } ix += incx; } } else { if (incy > 0) ky = 0; else ky = 1 - (n - 1) * incy - 1; ix = kx; for (i = 0; i < m; ++i) { if (x[ix] != 0.0) { temp = alpha * x[ix]; jy = ky; for (j = 0; j < n; ++j) { A(i, j) += temp * y[jy]; jy += incy; } } ix += incx; } } } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f06pact.c * * Copyright 1994 Numerical Algorithms Group. * * NAG C Library. * * dgemv level 2 C Blas. * * Author : George Levy, NAG Ltd., Oxford, U.K. * * Mark 3 re-issue, 1994. * */ #ifdef NAG_PROTO void f06pac(MatrixTranspose Trans, long m, long n, double alpha, const double a[], long tda, const double x[], long incx, double beta, double y[], long incy) #else void f06pac(Trans, m, n, alpha, a, tda, x, incx, beta, y, incy) MatrixTranspose Trans; long m, n; double alpha, a[]; long tda; double x[]; long incx; double beta, y[]; long incy; #endif { /* Local variables */ long info; double temp; long lenx, leny, i, j; long ix, iy, jx, jy, kx, ky; /* * Purpose * ======= * * dgemv performs one of the matrix-vector operations * * Y := ALPHA*A*X + BETA*Y, or Y := ALPHA*A'*X + BETA*Y, * * where ALPHA and BETA are scalars, X and Y are vectors and A is an * m by n matrix. * * Parameters * ========== * * Trans MatrixTranspose - specifies the operation to be performed as * follows: * Trans == NoTranspose Y := ALPHA*A*X + BETA*Y. * Trans == Transpose Y := ALPHA*A'*X + BETA*Y. * Trans == ConjugateTranspose Y := ALPHA*A'*X + BETA*Y. * * * m long - specifies the number of rows of the matrix A. * (m must be at least zero). * * * n long - specifies the number of columns of the matrix A. * (n must be at least zero). * * * alpha double - specifies the scalar ALPHA. * * * a const double [] - a[0] contains A(0, 0), the element in the first row * and first column of matrix A. The mapping used from a to A is * A(I, J) = a[(J) + (I) * tda] where I is the row index, with range * I = 0 .. m - 1, and J is the column index, with range J = 0 .. n - 1. * * * tda long - specifies the trailing dimension used in the definition * of matrix A. (tda must be at least MAX( 1, n)). * * * x const double [] - x is a pointer to a vector of dimension at * least ( 1 + ( n - 1 )*ABS( incx ) ) when Trans = NoTranspose * and at least ( 1 + ( m - 1 )*ABS( incx ) ) otherwise. * The first element of the vector is contained in x[0], and the (i+1)th * in x[i]. The incremented array x must contain the vector X. * * * incx long - specifies the increment for the elements of * X. (incx must not be zero). * * * beta double - specifies the scalar BETA. When beta is * supplied as zero then y need not be set on input. * * * y double [] - y is a pointer to a vector of dimension at least * ( 1 + ( m - 1 )*ABS( incy ) ) when Trans = NoTranspose * and at least ( 1 + ( n - 1 )*ABS( incy ) ) otherwise. * The first element of the vector is contained in y[0], and the (i+1)th * in y[i]. With beta non-zero, the incremented array y must contain * the vector Y. * On exit, y is overwritten by the updated vector y. * * * incy long - specifies the increment for the elements of the vector * Y. (incy must not be zero). * */ #define A9(I,J) a[(J)+(I) * tda] info = 0; if (Trans != NoTranspose && Trans != Transpose && Trans != ConjugateTranspose) info = 1; else if (m < 0) info = 2; else if (n < 0) info = 3; else if (tda < MAX(1,n)) info = 6; else if (incx == 0) info = 8; else if (incy == 0) info = 11; if (info != 0) { f06aaz("f06pac/dgemv", info); return; } /* Quick return if possible. */ if (m == 0 || n == 0 || (alpha == 0.0 && beta == 1.0)) return; /* * Set lenx and leny, the lengths of the vectors X and Y, and set * up the start points in x and y. */ if (Trans == NoTranspose) { lenx = n; leny = m; } else { lenx = m; leny = n; } if (incx > 0) kx = 0; else kx = - (lenx - 1) * incx; if (incy > 0) ky = 0; else ky = - (leny - 1) * incy; /* * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form Y := BETA*Y. */ if (beta != 1.0) { if (incy == 1) { if (beta == 0.0) for (i = 0; i < leny; ++i) y[i] = 0.0; else { for (i = 0; i < leny; ++i) y[i] *= beta; } } else { iy = ky; if (beta == 0.0) { for (i = 0; i < leny ; ++i) { y[iy] = 0.0; iy += incy; } } else { for (i = 0; i < leny; ++i) { y[iy] *= beta; iy += incy; } } } } if (alpha == 0.0) return; if (Trans==NoTranspose) { /* Form Y := ALPHA*A*X + Y. */ iy = ky; if (incx == 1) { for (i = 0; i < m; ++i) { temp = 0.0; for (j = 0; j < n; ++j) { temp += x[j] * A9(i,j); } y[iy] += temp * alpha; iy += incy; } } else { for (i = 0; i < m; ++i) { jx = kx; temp = 0.0; for (j = 0; j < n; ++j) { temp += x[jx] * A9(i,j); jx += incx; } y[iy] += temp * alpha; iy += incy; } } } else { /* Form Y := ALPHA*A'*X + Y. */ ix = kx; if (incy == 1) { for (i = 0; i < m; ++i) { if (x[ix] != 0.0) { temp = alpha * x[ix]; for (j = 0; j < n; ++j) { y[j] += temp * A9(i,j); } } ix += incx; } } else { for (i = 0; i < m; ++i) { if (x[ix] != 0.0) { temp = alpha * x[ix]; jy = ky; for (j = 0; j < n; ++j) { y[jy] += temp * A9(i,j); jy += incy; } } ix += incx; } } } } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f06pcct.c * * Copyright 1994 Numerical Algorithms Group. * * NAG C Library. * * dsymv level 2 C Blas. * * Author : George Levy, NAG Ltd., Oxford, U.K. * * Mark 3 re-issue, 1994. * */ #ifdef NAG_PROTO void f06pcc(MatrixTriangle UpperLower, long n, double alpha, const double a[], long tda, const double x[], long incx, double beta, double y[], long incy) #else void f06pcc(UpperLower, n, alpha, a, tda, x, incx, beta, y, incy) MatrixTriangle UpperLower; long n; double alpha; double a[]; long tda; double x[]; long incx; double beta; double y[]; long incy; #endif { /* Local variables */ long info; double temp1, temp2; long i, j; long ix, iy, jx, jy, kx, ky; #define A(I,J) a[(J)+(I)*tda] /* Purpose * ======= * * dsymv performs the matrix-vector operation * * Y := ALPHA*A*X + BETA*Y, * * where ALPHA and BETA are scalars, X and Y are n element vectors * and A is an n by n symmetric matrix. * * Parameters * ========== * * UpperLower MatrixTriangle - specifies whether the upper or lower * triangular part of the matrix A is to be referenced as follows: * * UpperLower == UpperTriangle Only the upper triangular part of A * is to be referenced. * * UpperLower == LowerTriangle Only the lower triangular part of A * is to be referenced. * * * n long - specifies the order of the matrix A. * (n must be at least zero). * * * alpha double - specifies the scalar ALPHA. * * * a const double [] - a[0] contains A(0, 0), the element in the first * row and first column of matrix A. The mapping used from a to A is * A(I, J) = a[(J) + (I) * tda] where I is the row index, with range * I = 0 .. n - 1, and J is the column index, with range J = 0 .. n - 1. * * With UpperLower = UpperTriangle, the n by n upper triangular part of * the matrix A must contain the upper triangular part of the symmetric * matrix and the strictly lower triangular part of A is not referenced. * * With UpperLower = LowerTriangle, the n by n lower triangular part of * the matrix A must contain the lower triangular part of the symmetric * matrix and the strictly upper triangular part of A is not referenced. * * * tda long - specifies the trailing dimension used in the definition * of matrix A. (tda must be at least MAX( 1, n)). * * x const double [] - x is a pointer to a vector of dimension at least * ( 1 + ( n - 1 )*ABS(incx) ). The first element of the vector is contained * in x[0], and the (i+1)th in x[i]. The incremented array x must contain * the n element vector X. * * * incx long - specifies the increment for the elements of x. * (incx must not be zero). * * * beta double - specifies the scalar BETA. When beta is supplied as zero * then y need not be set on input. * * * y double [] - y is a pointer to a vector of dimension at least * ( 1 + ( n - 1 )*ABS(incy) ). The first element of the vector is contained * in y[0], and the (i+1)th in y[i]. * On entry, the incremented array y must contain the n element vector Y. * On exit, y[] is overwritten by the updated vector y. * * * incy long - specifies the increment for the elements of * y. (incy must not be zero). * * */ info = 0; if (UpperLower!=UpperTriangle && UpperLower!=LowerTriangle) info = 1; else if (n < 0) info = 2; else if (tda < MAX(1, n)) info = 5; else if (incx == 0) info = 7; else if (incy == 0) info = 10; if (info != 0) { f06aaz("f06pcc/dsymv", info); return; } /* Quick return if possible. */ if (n == 0 || (alpha == 0.0 && beta == 1.0)) return; /* Set up the start points in x and y. */ if (incx > 0) kx = 1; else kx = 1 - (n - 1) * incx; if (incy > 0) ky = 1; else ky = 1 - (n - 1) * incy; /* * Start the operations. In this version the elements of A are * accessed sequentially with one pass through the triangular part * of A. */ /* First form Y := BETA*Y. */ if (beta != 1.0) { if (incy == 1) if (beta == 0.0) for (i = 0; i < n; ++i) y[i] = 0.0; else for (i = 0; i < n; ++i) y[i] = beta * y[i]; else { iy = ky - 1; if (beta == 0.0) for (i = 0; i < n; ++i) { y[iy] = 0.0; iy += incy; } else for (i = 0; i < n; ++i) { y[iy] = beta * y[iy]; iy += incy; } } } if (alpha == 0.0) return; if (UpperLower==UpperTriangle) { /* Form Y when A is stored in upper triangle. */ if (incx == 1 && incy == 1) { for (i = 0; i < n; ++i) { temp1 = alpha * x[i]; temp2 = 0.0; for (j = i + 1; j < n; ++j) { y[j] += temp1 * A(i, j); temp2 += A(i, j) * x[j]; } y[i] += temp1 * A(i, i) + alpha * temp2; } } else { ix = kx - 1; iy = ky - 1; for (i = 0; i < n; ++i) { temp1 = alpha * x[ix]; temp2 = 0.0; jy = iy; jx = ix; for (j = i + 1; j < n; ++j) { jx += incx; jy += incy; y[jy] += temp1 * A(i, j); temp2 += A(i, j) * x[jx]; } y[iy] += temp1 * A(i, i) + alpha * temp2; iy += incy; ix += incx; } } } else { /* Form Y when A is stored in lower triangle. */ if (incx == 1 && incy == 1) { for (i = 0; i < n; ++i) { temp1 = alpha * x[i]; temp2 = 0.0; for (j = 0; j < i; ++j) { y[j] += temp1 * A(i, j); temp2 += A(i, j) * x[j]; } y[i] += temp1 * A(i,i) + alpha * temp2; } } else { ix = kx - 1; iy = ky - 1; for (i = 0; i < n; ++i) { temp1 = alpha * x[ix]; temp2 = 0.0; jy = ky - 1; jx = kx - 1; for (j = 0; j < i; ++j) { y[jy] += temp1 * A(i, j); temp2 += A(i, j) * x[jx]; jx += incx; jy += incy; } y[iy] += temp1 * A(i,i) + alpha * temp2; iy += incy; ix += incx; } } } } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f06aaz.c * * Mark 1 Release. Copyright 1989 Numerical Algorithms Group. * Mark 2 Global Revision 1991 * Nag/ removed from include file names. SPD. * * NAG C Library * * Derived from: * * MARK 12 RELEASE. NAG COPYRIGHT 1986. * * Purpose * ======= * * F06AAZ is an error handler for the Level 2 BLAS routines. * * It is called by the Level 2 BLAS routines if an input parameter is * invalid. * * Parameters * ========== * * SRNAME - CHARACTER*13. * On entry, SRNAME specifies the name of the routine which * called F06AAZ. * * INFO - long. * On entry, INFO specifies the position of the invalid * parameter in the parameter-list of the calling routine. * * * Auxiliary routine for Level 2 Blas. * * Written on 20-July-1986. * * f06aaz.f -- translated by f2c (version 0.6, all of Fortran 77). * * Revised Mark 2 to replace explicit initialisation of error structure * by the macro INIT_FAIL. SPD 19.7.91 * * Tranlated 16 August 1989 */ #ifdef NAG_PROTO void f06aaz(char *srname, long info) #else void f06aaz(srname, info) char *srname; long info; #endif { NagError fail; char buf[NAG_ERROR_BUF_LEN]; INIT_FAIL(fail); fail.print = TRUE; Vsprintf(buf, "NE_F06_BAD_PARAM:\n\ On entry parameter number %1ld had an illegal value.", info); printf("should be errorhandler here! %s\n",buf); /*p01acc(buf, NE_F06_BAD_PARAM, srname, &fail);*/ } double x02amc(void) { return DBL_MIN; /* x02amc = 2.225073858507201400e-308 */ } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f01ajct.c * * Copyright 1990 Numerical Algorithms Group Ltd. * * NAG C Library * * Reduces the given lower triangle of a * symmetric matrix, A, stored in the array A[n][n], to * tridiagonal form using Householders reduction. * * Derived from: NAG Fortran Library F01AJF(Last Revision Mark 14) * * C version: Shah Datardina, August 1990, NAG Ltd, Oxford, U.K. * * Mark 1, 1990 * Mark 2 Global Revision 1991 * Nag/ removed from include file names. SPD. * Array pointers made explicit. GOM 17.12.91 * Mark 3a revised. IER-1694. (Apr 1994). * A Mark 1 coding mistake rectified. * Fortran Mark15a revision incorporated. SPD 22/4/94. */ #ifdef NAG_PROTO void f01ajc(long n, double a[], long tda, double d[], double e[], double z[], long tdz) #else void f01ajc(n, a, tda, d, e, z, tdz) long n; double a[]; long tda; double d[], e[], z[]; long tdz; #endif { #define A8(I,J) a[(I)*tda + (J)] #define Z(I,J) z[(I)*tdz + (J)] /* Local variables */ double f, g, h; long i, j, k, l; double scale, small; double hh; long ii; small = X02AMC; for (i=0; i small) { f06edc(l,1.0/scale, d, (long)1); h = f06eac(l, d, (long)1, d, (long)1); f = d[i-2]; g = sqrt(h); if (f>=0.0) g = -g; e[i-1] = g*scale; h -= f*g; d[i-2] = f-g; /* Copy u */ for (j=0; j m || right && k2 > n) return; if (left) { if (pivot==VariablePivot) { if (direct==ForwardSequence) { for (j = 0; j < n; ++j) { aij = A10(k1-1,j); for (i = k1-1; i < k2-1; ++i) { temp = A10(i+1,j); A10(i,j) = s[i] * temp + c[i] * aij; aij = c[i] * temp - s[i] * aij; } A10(k2-1,j) = aij; } } else if (direct==BackwardSequence) { for (j = 0; j < n; ++j) { aij = A10(k2-1,j); for (i = k2-2; i >= k1-1; --i) { temp = A10(i,j); A10(i+1,j) = c[i] * aij - s[i] * temp; aij = s[i] * aij + c[i] * temp; } A10(k1-1,j) = aij; } } } else if (pivot==TopPivot) { if (direct==ForwardSequence) { for (j = 0; j < n; ++j) { temp = A10(k1-1,j); for (i = k1-1; i < k2-1; ++i) { aij = A10(i+1,j); A10(i+1,j) = c[i] * aij - s[i] * temp; temp = s[i] * aij + c[i] * temp; } A10(k1-1,j) = temp; } } else if (direct==BackwardSequence) { for (j = 0; j < n; ++j) { temp = A10(k1-1,j); for (i = k2-2; i >= k1-1; --i) { aij = A10(i+1,j); A10(i+1,j) = c[i] * aij - s[i] * temp; temp = s[i] * aij + c[i] * temp; } A10(k1-1,j) = temp; } } } else if (pivot==BottomPivot) { if (direct==ForwardSequence) { for (j = 0; j < n; ++j) { temp = A10(k2-1,j); for (i = k1; i < k2; ++i) { aij = A10(i,j); A10(i,j) = s[i] * temp + c[i] * aij; temp = c[i] * temp - s[i] * aij; } A10(k2-1,j) = temp; } } else if (direct==BackwardSequence) { for (j = 0; j < n; ++j) { temp = A10(k2-1,j); for (i = k2-2; i >= k1-1; --i) { aij = A10(i,j); A10(i,j) = s[i] * temp + c[i] * aij; temp = c[i] * temp - s[i] * aij; } A10(k2-1,j) = temp; } } } } else if (right) { if (pivot==VariablePivot) { if (direct==ForwardSequence) { for (j = k1-1; j < k2-1; ++j) { if (c[j] != 1.0 || s[j] != 0.0) { ctemp = c[j]; stemp = s[j]; for (i = 0; i < m; ++i) { temp = A10(i,j+1); A10(i,j+1) = ctemp * temp - stemp * A10(i,j); A10(i,j) = stemp * temp + ctemp * A10(i,j); } } } } else if (direct==BackwardSequence) { for (j = k2-2; j >= k1-1; --j) { if (c[j] != 1.0 || s[j] != 0.0) { ctemp = c[j]; stemp = s[j]; for (i = m-1; i >= 0; --i) { temp = A10(i,j+1); A10(i,j+1) = ctemp*temp - stemp*A10(i,j); A10(i,j) = stemp*temp + ctemp*A10(i,j); } } } } } else if (pivot==TopPivot) { if (direct==ForwardSequence) { for (j = k1; j < k2; ++j) { ctemp = c[j - 1]; stemp = s[j - 1]; if (ctemp != 1.0 || stemp != 0.0) { for (i = 0; i < m; ++i) { temp = A10(i,j); A10(i,j) = ctemp*temp - stemp*A10(i,k1-1); A10(i,k1-1) = stemp*temp + ctemp*A10(i,k1-1); } } } } else if (direct==BackwardSequence) { for (j = k2-1; j >= k1; --j) { ctemp = c[j - 1]; stemp = s[j - 1]; if (ctemp != 1.0 || stemp != 0.0) { for (i = m-1; i >= 0; --i) { temp = A10(i,j); A10(i,j) = ctemp*temp - stemp*A10(i,k1-1); A10(i,k1-1) = stemp*temp + ctemp*A10(i,k1-1); } } } } } else if (pivot==BottomPivot) { if (direct==ForwardSequence) { for (j = k1-1; j < k2-1; ++j) { if (c[j] != 1.0 || s[j] != 0.0) { ctemp = c[j]; stemp = s[j]; for (i = 0; i < m; ++i) { temp = A10(i,j); A10(i,j) = stemp*A10(i,k2-1) + ctemp*temp; A10(i,k2-1) = ctemp*A10(i,k2-1) - stemp*temp; } } } } else if (direct==BackwardSequence) { for (j = k2-2; j >= k1-1; --j) { if (c[j] != 1.0 || s[j] != 0.0) { ctemp = c[j]; stemp = s[j]; for (i = m-1; i >= 0; --i) { temp = A10(i,j); A10(i,j) = stemp*A10(i,k2-1) + ctemp*temp; A10(i,k2-1) = ctemp*A10(i,k2-1) - stemp*temp; } } } } } } } /* f06qxc */ /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f02amct.c * * Copyright 1990 Numerical Algorithms Group Ltd. * * NAG C Library * * Finds the eigenvalues and eigenvectors of a * tridiagonal matrix * * Derived from: NAG Fortran Library F02AMF(Last Revision Mark 14) * * C version: Shah Datardina, August 1990, NAG Ltd, Oxford, U.K. * * Mark 1, 1990 * Mark 2 Revised 1991 * fail.message now set to "NE_NOERROR" when no error found. SPD * Mark 2 Global Revision 1991 * Nag/ removed from include file names. SPD. * Mark 2 Revised. (GOM - January 1991). * character variable type replaced by Enum variable types. * */ #ifdef NAG_PROTO void f02amc(long n, double eps, double *d, double *e, double *z, long tdz, NagError *fail) #else void f02amc(n, eps, d, e, z, tdz, fail) long n; double eps; double *d, *e, *z; long tdz; NagError *fail; #endif { #define Z(I,J) z[(I)*tdz + (J)] #define VLEN 128 /* Local variables */ long iseg, ipos; double b, c, f, g, h; long i, j, k, l, m; double p; double r, s; long i1; double cc[VLEN], ss[VLEN]; char buf[NAG_ERROR_BUF_LEN]; if (fail) fail->code = NE_NOERROR; if (n!=1) { for (i=1; i b; ++m) ; if (m != l+1) { do { if (j<=0) { Vsprintf(buf,"NE_TOO_MANY_ITERATIONS:\n\ More than %1ld iterations are required to isolate\ all the eigenvalues",30*n); printf("In f02amc: %s\n",buf); /*p01acc(buf, NE_TOO_MANY_ITERATIONS,"f02amc",fail);*/ return; } --j; /* Form shift */ g = d[l]; h = d[l + 1] - g; if (FABS(h) < FABS(e[l])) { p = h * 0.5 / e[l]; r = sqrt(p * p + 1.0); h = p + r; if (p<0.0) { h = p - r; } d[l] = e[l] / h; } else { p = e[l] * 2.0 / h; r = sqrt(p * p + 1.0); d[l] = e[l] * p / (r + 1.0); } h = g - d[l]; i1 = l+1 + 1; if (i1<=n) { for (i=i1; i<=n; ++i) { d[i-1] -= h; } } f += h; /* QL transformation */ p = d[m-1]; c = 1.0; s = 0.0; for (k = m-1; k >= l+1; k += -VLEN) { iseg = MAX(k-VLEN+1,l+1); for (i=k; i>=iseg; --i) { g = c * e[i-1]; h = c * p; if (FABS(p) >= FABS(e[i-1])) { c = e[i-1] / p; r = sqrt(c * c + 1.0); e[i] = s * p * r; s = c / r; c = 1.0 / r; } else { c = p / e[i-1]; r = sqrt(c * c + 1.0); e[i] = s * e[i-1] * r; s = 1.0 / r; c /= r; } p = c * d[i-1] - s * g; d[i] = h + s * (c * g + s * d[i-1]); /* Store rotations */ cc[VLEN - k + i - 1] = c; ss[VLEN - k + i - 1] = -s; } /* Update vectors */ ipos = VLEN - k + iseg; f06qxc(RightSide, VariablePivot, BackwardSequence, n, (long)(k-iseg+2), (long)1, (long)(k-iseg+2), &cc[ipos - 1], &ss[ipos - 1], &Z(0,iseg - 1), tdz); } e[l] = s * p; d[l] = c * p; } while (FABS(e[l]) > b) ; } d[l] += f; } /* Order eigenvalues and eigenvectors */ for (i=0; icode==NE_NOERROR) printf("NE_NOERROR: No error\n");*/ /* Library Wide Error Messges */ /*p01acc((char *)0, NE_NOERROR, "f02amc", fail);*/ return; } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /*------------------------------------------------------------------*/ void f02abc(long n, double *a, long tda, double *r, double *v, long tdv, NagError *fail) /* Eigenvalues and eigenvectors of a real symmetrix matrix. */ /*------------------------------------------------------------------*/ { double *e; /* Work Array */ char buf[NAG_ERROR_BUF_LEN]; NagError savefail; INIT_FAIL(savefail); if (fail) fail->code = NE_NOERROR; if (n < 1) { Vsprintf(buf,"NE_INT_ARG_LT:\n\ On entry, %s must not be less than %s: %s = %ld.", "n", "1", "n", n); printf("In %s: %s\n","f02abc", buf); /*p01acc(buf, NE_INT_ARG_LT, "f02abc", fail);*/ return; } if (tda < n) { Vsprintf(buf,"NE_2_INT_ARG_LT:\n\ On entry %s = %2ld while %s = %2ld.\n\ These parameters must satisfy %s >= %s", "tda", tda, "n", n, "tda", "n"); printf("In %s: %s\n","f02abc", buf); /*p01acc(buf, NE_2_INT_ARG_LT, "f02abc", fail);*/ return; } if (tdv < n) { Vsprintf(buf,"NE_2_INT_ARG_LT:\n\ On entry %s = %2ld while %s = %2ld.\n\ These parameters must satisfy %s >= %s", "tdv", tdv, "n", n, "tdv", "n"); printf("In %s: %s\n","f02abc", buf); /*p01acc(buf, NE_2_INT_ARG_LT, "f02abc", fail);*/ return; } e = NAG_ALLOC(n, double); if (e==(double *)0) { printf("In %s: %s\n","f02abc", "NE_ALLOC_FAIL:\n\ Memory allocation failed."); /*p01acc((char *)0, NE_ALLOC_FAIL, "f02abc", fail);*/ return; } f01ajc(n, a, tda, r, e, v, tdv); f02amc(n, X02AJC, r, e, v, tdv, &savefail); if (savefail.code != NE_NOERROR) { NAG_FREE(e); Vsprintf(buf,"NE_TOO_MANY_ITERATIONS:\n\ More than %1ld iterations are required to isolate all the eigenvalues", 30*n); printf("In %s: %s\n","f02abc", buf); /*p01acc(buf, NE_TOO_MANY_ITERATIONS,"f02abc",fail);*/ } NAG_FREE(e); /* if (fail && fail->code==NE_NOERROR) printf("In %s: %s\n","f02abc","NE_NOERROR: No error");*/ /*p01acc((char *)0, NE_NOERROR, "f02abc", fail);*/ } /*------------------------------------------------------------------*/ /* New Routines for ASYMMETRIC diagonalization */ /*------------------------------------------------------------------*/ /* a02bact.c * * Copyright 1991 Numerical Algorithms Group * * NAG C Library * * Form complex number from real and imaginary parts * * Mick Pont, NAG Ltd, Oxford, U.K., July 1991. * * Mark 2, 1991 */ #ifdef NAG_PROTO Complex a02bac(double x, double y) #else Complex a02bac(x, y) double x, y; #endif { Complex z; z.re = x; z.im = y; return z; } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* a02bcct.c * * Copyright 1991 Numerical Algorithms Group * * NAG C Library * * Take imaginary part of complex number * * Mick Pont, NAG Ltd, Oxford, U.K., July 1991. * * Mark 2, 1991 */ #ifdef NAG_PROTO double a02bcc(Complex z) #else double a02bcc(z) Complex z; #endif { return z.im; } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* a02bbct.c * * Copyright 1991 Numerical Algorithms Group * * NAG C Library * * Take real part of complex number * * Mick Pont, NAG Ltd, Oxford, U.K., July 1991. * * Mark 2, 1991 */ #ifdef NAG_PROTO double a02bbc(Complex z) #else double a02bbc(z) Complex z; #endif { return z.re; } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* a02cdct.c * * Copyright 1991 Numerical Algorithms Group * * NAG C Library * * Complex division * * Derived from: NAG Fortran Library A02ACF (Last Revision Mark 11.5) * * C version: Ian Hounam, NAG Ltd, Oxford, U.K., September 1989. * * Renamed a02cdc, and parameter names changed: Mick Pont, July 1991. * * Mark 2, 1991 */ #ifdef NAG_PROTO Complex a02cdc(Complex z1, Complex z2) #else Complex a02cdc(z1, z2) Complex z1, z2; #endif { double a, h; Complex z; if (FABS(z2.re) <= FABS(z2.im)) { h = z2.re / z2.im; a = 1.0 / (h * z2.re + z2.im); z.re = (h * z1.re + z1.im) * a; z.im = (h * z1.im - z1.re) * a; } else { h = z2.im / z2.re; a = 1.0 / (h * z2.im + z2.re); z.re = (z1.re + (h * z1.im)) * a; z.im = (z1.im - (h * z1.re)) * a; } return z; } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* a02dbct.c * * Copyright 1991 Numerical Algorithms Group * * NAG C Library * * Returns the absolute value of a complex number * * Derived from: NAG Fortran Library A02ABF (Last Revision Mark 11.5) * * C version: Ian Hounam, NAG Ltd, Oxford, U.K., September 1989. * * Renamed a02dbc, and parameter name changed: Mick Pont, July 1991. * * Mark 2, 1991 */ #ifdef NAG_PROTO double a02dbc(Complex z) #else double a02dbc(z) Complex z; #endif { double h, zi, zr, pwr; zr = FABS(z.re); zi = FABS(z.im); if (zi > zr) { h = zr; zr = zi; zi = h; } if (zi != 0.0) { pwr = zi / zr; h = zr * sqrt(1.0 + pwr * pwr); return h; } else return zr; } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* x02akc.c * * Mark 2 Release. * Copyright 1991 Numerical Algorithms Group Ltd., Oxford, U.K. * * This function was generated by getcons. * */ double x02akc(void) { return DBL_MIN; /* x02akc = 2.225073858507201400e-308 */ } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f06pjct.c * * Copyright 1994 Numerical Algorithms Group. * * NAG C Library. * * dtrsv level 2 C Blas. * * Author : George Levy, NAG Ltd., Oxford, U.K. * * Mark 3 re-issue, 1994. * */ #ifdef NAG_PROTO void f06pjc(MatrixTriangle UpperLower, MatrixTranspose Trans, MatrixUnitTriangular TriangularMatrix, long n, const double a[], long tda, double x[], long incx) #else void f06pjc(UpperLower, Trans, TriangularMatrix, n, a, tda, x, incx) MatrixTriangle UpperLower; MatrixTranspose Trans; MatrixUnitTriangular TriangularMatrix; long n; double a[]; long tda; double x[]; long incx; #endif { /* Local variables */ long info; double temp; long i, j; long ix, jx, kx; Boolean nounit; /* * Purpose * ======= * * dtrsv solves one of the systems of equations * * A*X = B, or A'*X = B, * * where B and X are n element vectors and A is an n by n unit, or * non-unit, upper or lower triangular matrix. * * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. * * Parameters * ========== * * UpperLower MatrixTriangle - specifies whether the matrix A is an * upper or lower triangular matrix. * * * Trans MatrixTranspose - specifies the equations to be solved as * follows: * Trans = NoTranspose A*X = B. * Trans = Transpose A'*X = B. * Trans = ConjugateTranspose A'*X = B. * * * TriangularMatrix MatrixUnitTriangular - specifies whether or not A is * unit triangular as follows: * TriangularMatrix = UnitTriangular A is assumed to be unit triangular. * TriangularMatrix = NotUnitTriangular A is not assumed to be unit * triangular. * * * n long - specifies the order of the matrix A. * (n must be at least zero). * * * a const double [] - a[0] contains A(0, 0) the element in the first row * and first column of matrix A. The mapping used from a to A is * A(I, J) = a[(J) + (I) * tda] where I is the row index and J is the column * index. * With UpperLower = UpperTriangular, the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * With UpperLower = LowerTriangular, the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when TriangularMatrix = UnitTriangular, the diagonal elements of * A are not referenced either, but are assumed to be unity. * * * tda long - specifies the trailing dimension used in the definition * of matrix A. (tda must be at least MAX(1, n)). * * * x double [] - x is a pointer to a vector of dimension at least * (1+(n-1)*ABS(incx)). The first element of the vector is contained in x[0], * and the (i+1)th in x[i]. Before entry, the incremented array x must * contain the n element right-hand side vector B. On exit, x is overwritten * with the solution vector x. * * * incx long - incx specifies the increment for the elements of * x. (incx must not be zero). * * */ #define A(I,J) a[(J)+(I)*tda] info = 0; if (UpperLower!=UpperTriangle && UpperLower!=LowerTriangle) info = 1; else if (Trans!=NoTranspose && Trans!=Transpose && Trans!=ConjugateTranspose) info = 2; else if (TriangularMatrix!=UnitTriangular && TriangularMatrix!=NotUnitTriangular) info = 3; else if (n<0) info = 4; else if (tda= 0; --i) { temp = x[i]; for (j = i + 1; j < n ; ++j) temp -= A(i, j) * x[j]; if (nounit) temp /= A(i, i); x[i] = temp; } } else { ix = kx + (n - 1) * incx - 1; for (i = n - 1; i >= 0; --i) { temp = x[ix]; jx = ix; for (j = i + 1; j < n ; ++j) { jx += incx; temp -= A(i, j) * x[jx]; } if (nounit) temp /= A(i, i); x[ix] = temp; ix -= incx; } } } else { if (incx == 1) { for (i = 0; i < n; ++i) { temp = x[i]; for (j = 0; j < i; ++j) temp -= A(i, j) * x[j]; if (nounit) temp /= A(i, i); x[i] = temp; } } else { ix = kx - 1; for (i = 0; i < n; ++i) { temp = x[ix]; jx = kx - 1; for (j = 0; j < i; ++j) { temp -= A(i, j) * x[jx]; jx += incx; } if (nounit) temp /= A(i, i); x[ix] = temp; ix += incx; } } } } else { /* Form X := inv(A')*X. */ if (UpperLower==UpperTriangle) { if (incx == 1) { for (i = 0; i < n; ++i) { if (nounit) x[i] /= A(i, i); if (x[i] != 0.0) { temp = x[i]; for (j = i + 1; j < n; ++j) x[j] -= A(i, j) * temp; } } } else { ix = kx - 1; for (i = 0; i < n; ++i) { if (nounit) x[ix] /= A(i, i); jx = ix; if (x[ix] != 0.0) { temp = x[ix]; for (j = i + 1; j < n; ++j) { jx += incx; x[jx] -= A(i, j) * temp; } } ix += incx; } } } else { if (incx == 1) { for (i = n - 1; i >= 0; --i) { if (nounit) x[i] /= A(i, i); if (x[i] != 0.0) { temp = x[i]; for (j = 0; j < i; ++j) x[j] -= A(i, j) * temp; } } } else { ix = kx + (n - 1) * incx - 1; for (i = n - 1; i >= 0; --i) { if (nounit) x[ix] /= A(i, i); jx = kx - 1; if (x[ix] != 0.0) { temp = x[ix]; for (j = 0; j < i; ++j) { x[jx] -= A(i, j) * temp; jx += incx; } } ix -= incx; } } } } } /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f01atzt.c * * Copyright 1990 Numerical Algorithms Group * * NAG C Library * * Interchanges elements 1 to k of columns j and m, * and elements l to n of rows j and m of matrix A. * * Derived from: NAG Fortran Library F01ATZ (Last Revision Mark 11.5) * * C version: Shah Datardina, August 1990, NAG Ltd, Oxford, U.K. * * Mark 1, 1990 * Mark 2 Global Revision 1991 * Nag/ removed from include file names. SPD. */ #ifdef NAG_PROTO void f01atz(long m, double *a, long tda, double *d, long k, long l, long n, long j) #else void f01atz(m, a, tda, d, k, l, n, j) long m; double *a; long tda; double *d; long k, l; long n; long j; #endif { #define A1(I,J) a[(I)*tda + (J)] /* Local variables */ double f; long i; d[m-1] = (double) (j); if (j!=m) { for (i=0; i k) { goto L180; } /* if */ for (j = l; j <= k; ++j) { c = 0.0; for (i = l; i <= k; ++i) { if (i == j) { goto L120; } /* if */ c += FABS(A2(i, j)); L120: ;} /* for */ if (c == 0.0) { goto L160; } /* if */ } /* for */ goto L180; L160: f01atz(l, a, tda, d, k, l, n, j); ++l; goto L100; /* Now balance the submatrix in rows l through k. */ L180: *low = l; *lhi = k; if (l > k) { goto L220; } /* if */ for (i = l; i <= k; ++i) { D(i) = 1.0; } /* for */ L220: noconv = FALSE; if (l > k) { goto L420; } /* if */ for (i = l; i <= k; ++i) { c = 0.0; r = 0.0; for (j = l; j <= k; ++j) { if (j == i) { goto L240; } /* if */ c += FABS(A2(j, i)); r += FABS(A2(i, j)); L240: ;} /* for */ g = r / (double) (ib); f = 1.0; s = c + r; L260: if (c >= g) { goto L280; } /* if */ f *= (double) (ib); c *= b2; goto L260; L280: g = r * (double) (ib); L300: if (c < g) { goto L320; } /* if */ f /= (double) (ib); c /= b2; goto L300; L320: if ((c + r) / f >= s * 0.95) { goto L400; } /* if */ g = 1.0 / f; D(i) *= f; noconv = TRUE; if (l > n) { goto L360; } /* if */ for (j = l; j <= n; ++j) { A2(i, j) *= g; } /* for */ L360: for (j = 1; j <= k; ++j) { A2(j, i) *= f; } /* for */ L400: ;} /* for */ L420: if (noconv) { goto L220; } /* if */ return ; } /* f01atc */ /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f01akct.c * * Copyright 1994 Numerical Algorithms Group. * * NAG C Library * * Given the unsymmetric matrix, A, STORED IN THE ARRAY a[n*n], * this subroutine reduces the sub-matrix of order l - k + 1, * which starts at the element A(K,K) and finishes at the * element A(l,l), to hessenberg form, H, by the direct * method(AN = NH). The matrix H is overwritten on a with * details of the transformations (n) stored in the remaining * triangle under h and in elements k to l of the array * INTGER(N). * * Derived from: NAG Fortran Library F01AKF (Last Revision Mark 12) * * C version: Shah Datardina, NAG Ltd., Oxford, U.K. * * Mark 4 re-issue, 1994. * */ /* Table of constant values */ #ifdef NAG_PROTO void f01akc(long n, long k, long l, double a[], long tda, long intger[]) #else void f01akc(n, k, l, a, tda, intger) long n, k, l; double a[]; long tda, intger[]; #endif { /* Local variables */ double c_b10 = 1.0; double c_b19 = -1.0; double d__1; double x, y; long i__2, i__3; long i, j, m; long k1; /* Function Body */ #define INTGER(I) intger[(I)-1] #define A3(I,J) a[((I)-1) * ( tda) + ((J)-1)] k1 = k + 1; if (k1 > l) { return; } for (j = k1; j <= n; ++j) { m = j; x = 0.0; if (j > l) { goto L120; } for (i = j; i <= l; ++i) { if ((d__1 = A3(i,j-1), ABS(d__1)) <= ABS(x)) { goto L20; } x = A3(i,j-1); m = i; L20: ; } INTGER(j) = m; if (m == j) { goto L80; } /* INTERCHANGE ROWS AND COLUMNS OF A. */ for (i = k; i <= n; ++i) { y = A3(m,i); A3(m,i) = A3(j,i); A3(j,i) = y; /* L40: */ } for (i = 1; i <= l; ++i) { y = A3(i,m); A3(i,m) = A3(i,j); A3(i,j) = y; /* L60: */ } L80: if (x != 0.0 && j < l) { for (i = j + 1; i <= l; ++i) { A3(i,j-1) /= x; /* L100: */ } i__2 = l - j; f06pac(NoTranspose, l, i__2, c_b10, &A3(1,j+1), tda, &A3(j+1,j-1), tda, c_b10, &A3(1,j), tda); } L120: i__2 = j - k; f06pjc(LowerTriangle, NoTranspose, UnitTriangular, i__2, &A3(k+1,k), tda, &A3(k+1,j), tda); if (j < l) { i__2 = l - j; i__3 = j - k; f06pac(NoTranspose, i__2, i__3, c_b19, &A3(j+1,k), tda, &A3(k+1,j), tda, c_b10, &A3(j+1,j), tda); } /* L140: */ } return; } /* f01akc_ */ /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f01apct.c * * Copyright 1994 Numerical Algorithms Group. * * NAG C Library * * Derived from: NAG Fortran Library F01APF (Last Revision Mark 11.5) * * C version: Shah Datardina, NAG Ltd., Oxford, U.K. * * FORM the matrix of accumulated transformations in the array * v[n*n] from the information left by subroutine f01akc * below the upper hessenberg matrix, H, in the array h[n*n] * and in the integer array intger[n]. * * Mark 4 re-issue, 1994. * */ #ifdef NAG_PROTO void f01apc(long n, long low, long iupp, long intger[], double h[], long tdh, double v[], long tdv) #else void f01apc(n, low, iupp, intger, h, tdh, v, tdv) long n, low, iupp, intger[]; double h[]; long tdh; double v[]; long tdv; #endif { /* Local variables */ double x; long i, j, m; long i1, ii, low1; /* Function Body */ #define INTGER(I) intger[(I)-1] #define V(I,J) v[((I)-1) * ( tdv) + ((J)-1)] #define H(I,J) h[((I)-1) * ( tdh) + ((J)-1)] for (i = 1; i <= n; ++i) { for (j = 1; j <= n; ++j) { V(i,j) = 0.0; /* L20: */ } V(i,i) = 1.0; /* L40: */ } low1 = low + 1; if (low1 > iupp) { return; } for (ii = low1; ii <= iupp; ++ii) { i = low1 + iupp - ii; i1 = i - 1; if (low1 > i1) { goto L80; } for (j = low1; j <= i1; ++j) { V(i,j) = H(i,j-1); /* L60: */ } L80: m = INTGER(i); if (m == i) { goto L120; } for (j = low1; j <= iupp; ++j) { x = V(m,j); V(m,j) = V(i,j); V(i,j) = x; /* L100: */ } L120: ; } return; } /* f01apc */ /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f01auct.c * * Copyright 1994 Numerical Algorithms Group. * * NAG C Library * * Derived from: NAG Fortran Library F01AUF (Last Revision Mark 11.5) * * C version: Shah Datardina, NAG Ltd., Oxford, U.K. * * Mark 4 re-issue, 1994. * */ #ifdef NAG_PROTO void f01auc(long n, long low, long lhi, long m, double d[], double z[], long tdz) #else void f01auc(n, low, lhi, m, d, z, tdz) long n, low, lhi, m; double d[], z[]; long tdz; #endif { /* Local variables */ double s; long i, j, k; long ii, lhi1, low1; /* Function Body */ #define D(I) d[(I)-1] #define Z1(I,J) z[((I)-1) * ( tdz) + ((J)-1)] if (low > lhi) { goto L60; } for (i = low; i <= lhi; ++i) { s = D(i); /* Left-hand eigenvectors are back transformed if the * foregoing statement is replaced by s=1/D(I) */ for (j = 1; j <= m; ++j) { Z1(i,j) *= s; /* L20: */ } /* L40: */ } L60: i = low; low1 = low - 1; if (low1 < 1) { goto L120; } for (ii = 1; ii <= low1; ++ii) { --i; k = (long) D(i); if (k == i) { goto L100; } for (j = 1; j <= m; ++j) { s = Z1(i,j); Z1(i,j) = Z1(k,j); Z1(k,j) = s; /* L80: */ } L100: ; } L120: lhi1 = lhi + 1; if (lhi1 > n) { return; } for (i = lhi1; i <= n; ++i) { k = (long) D(i); if (k == i) { goto L160; } for (j = 1; j <= m; ++j) { s = Z1(i,j); Z1(i,j) = Z1(k,j); Z1(k,j) = s; /* L140: */ } L160: ; } return; } /* f01auc */ /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f02aqct.c * * Copyright 1994 Numerical Algorithms Group * * NAG C Library * * Finds the eigenvalues and eigenvectors of a real matrix * which has been reduced to upper Hessenberg form. * * Derived from: NAG Fortran Library F02AQF(Last Revision Mark 13 ) * * Mark 3 re-issue, 1994. * * Mark 3 C version rederived from Fortran to re-insert Fortran control * structures. SPD August 1993. * */ #ifdef NAG_PROTO void f02aqc(long n, long low, long upp, double machep, double *h, long tdh, double *vecs, long tdvecs, Complex *eigvals, long *cnt, NagError *fail) #else void f02aqc(n, low, upp, machep, h, tdh, vecs, tdvecs, eigvals, cnt, fail) long n; long low, upp; double machep; double *h; long tdh; double *vecs; long tdvecs; Complex *eigvals; long *cnt; NagError *fail; #endif { /* Local variables */ double norm; long i, j, k, l, m; double p, q, r, s, t, u, w, x, y, z; long i1, m2, m3, na, ii; double ra, sa; long en, jj, kk, ll, mm; double vi, vr; long na1; Boolean notlas; long en2, nhs, itn, its, low1, upp1; Complex tmpc1, tmpc2, tmpc3; char buf[NAG_ERROR_BUF_LEN]; /* Parameter adjustments */ #define H1(I,J) h[(I-1)*tdh + (J-1)] #define VECS(I,J) vecs[(I-1)*tdvecs + (J-1)] #define EIGVALS(I) eigvals[(I-1)] #define CNT(I) cnt[(I-1)] /* * * HQR2 * Finds the eigenvalues and eigenvectors of a real matrix * which has been reduced to upper hessenberg form in the array * h with the accumulated transformations stored in * the array vecs. The real and imaginary parts of the * eigenvalues are formed in the array eigvals and the * eigenvectors are formed in the array vecs where * only one complex vector, corresponding to the root with * positive imaginary part, is formed for a complex pair. low * and upp are two integers produced in balancing where * eigenvalues are isolated in positions 1 to low-1 and upp+1 * to n. If balancing is not used low=1, upp=n. machep is the * relative machine precision. The subroutine will fail if * all eigenvalues take more than 30*n iterations. * */ if (fail) fail->code = NE_NOERROR; /* Compute matrix norm */ norm = 0.0; k = 1; for (i = 1; i <= n; ++i) { for (j = k; j <= n; ++j) { norm += FABS(H1(i, j)); } /* for */ k = i; } /* for */ nhs = n * (n + 1) / 2 + n - 1; /* Isolated roots */ if (low <= 1) { goto L80; } /* if */ j = low - 1; for (i = 1; i <= j; ++i) { EIGVALS(i).re = H1(i, i); EIGVALS(i).im = 0.0; CNT(i) = 0; } /* for */ L80: if (upp >= n) { goto L120; } /* if */ j = upp + 1; for (i = j; i <= n; ++i) { EIGVALS(i).re = H1(i,i); EIGVALS(i).im = 0.0; CNT(i) = 0; } /* for */ L120: en = upp; t = 0.0; itn = n * 30; L140: if (en < low) { goto L880; } /* if */ its = 0; na = en - 1; /* Look for single small sub-diagonal element */ L160: if (low + 1 > en) { goto L200; } /* if */ low1 = low + 1; for (ll = low1; ll <= en; ++ll) { l = en + low1 - ll; s = FABS(H1(l - 1,l - 1)) + FABS(H1(l, l)); if (s < X02AKC / X02AJC) { s = norm / (double) nhs; } /* if */ if (FABS(H1(l, l - 1)) <= machep * s) { goto L220; } /* if */ } /* for */ L200: l = low; L220: x = H1(en,en); if (l == en) { goto L740; } /* if */ y = H1(na,na); w = H1(en, na) * H1(na, en); if (l == na) { goto L760; } /* if */ if (itn <= 0) { goto L1500; } /* if */ /* Form shift */ if (its != 10 && its != 20) { goto L280; } /* if */ t += x; if (low > en) { goto L260; } /* if */ for (i = low; i <= en; ++i) { H1(i, i) -= x; } /* for */ L260: s = FABS(H1(en, na)) + FABS(H1(na,en-2)) ; x = s * 0.75; y = x; w = s * s * -0.4375; L280: ++its; --itn; /* Look for two consecutive small sub-diagonal elements */ if (l > en - 2) { goto L320; } /* if */ en2 = en - 2; for (mm = l; mm <= en2; ++mm) { m = l + en2 - mm; z = H1(m, m); r = x - z; s = y - z; p = (r * s - w) / H1(m+1, m) + H1(m, m+1); q = H1(m+1,m+1) - z - r - s; r = H1(m+2,m+1); s = (FABS(p)) + (FABS(q)) + (FABS(r)); p /= s; q /= s; r /= s; if (m == l) { goto L320; } /* if */ if (FABS(H1(m, m-1)) * ((FABS(q)) + (FABS(r))) <= machep * (FABS(p)) * (FABS(H1(m - 1, m - 1)) + (FABS(z)) + FABS(H1(m + 1, m + 1)))) { goto L320; } /* if */ } /* for */ L320: m2 = m + 2; if (m2 > en) { goto L360; } /* if */ for (i = m2; i <= en; ++i) { H1(i, i - 2) = 0.0; } /* for */ L360: m3 = m + 3; if (m3 > en) { goto L400; } /* if */ for (i = m3; i <= en; ++i) { H1(i, i - 3) = 0.0; } /* for */ L400: if (m > na) { goto L720; } /* if */ for (k = m; k <= na; ++k) { notlas = k != na; if (k == m) { goto L420; } /* if */ p = H1(k, k - 1); q = H1(k + 1, k - 1); r = 0.0; if (notlas) { r = H1(k + 2, k - 1); } /* if */ x = (FABS(p)) + (FABS(q)) + (FABS(r)); if (x == 0.0) { goto L700; } /* if */ p /= x; q /= x; r /= x; L420: s = sqrt(p * p + q * q + r * r); if (p < 0.0) { s = -s; } /* if */ if (k != m) { goto L440; } /* if */ if (l != m) { H1(k, k-1) = -H1(k, k - 1); } /* if */ goto L460; L440: H1(k, k - 1 ) = -s * x; L460: p += s; x = p / s; y = q / s; z = r / s; q /= p; r /= p; /* Row modification */ if (notlas) { goto L500; } /* if */ for (j = k; j <= n; ++j) { p = H1(k, j) + q * H1(k+1, j); H1(k + 1, j) -= p * y; H1(k, j) -= p * x; } /* for */ goto L540; L500: for (j = k; j <= n; ++j) { p = H1(k, j) + q * H1((k + 1), j ) + r * H1(k + 2 ,j); H1(k + 2 ,j) -= p * z; H1(k + 1 ,j) -= p * y; H1(k, j) -= p * x; } /* for */ L540: j = en; if (k + 3 < en) { j = k + 3; } /* if */ /* Column modification */ if (notlas) { goto L580; } /* if */ for (i = 1; i <= j; ++i) { p = x * H1(i, k) + y * H1(i, k + 1); H1(i, k + 1 ) -= p * q; H1(i, k) -= p; } /* for */ goto L620; L580: for (i = 1; i <= j; ++i) { p = x * H1(i, k) + y * H1(i , k + 1) + z * H1(i , k + 2); H1(i , k + 2 ) -= p * r; H1(i , k + 1) -= p * q; H1(i, k) -= p; } /* for */ /* Accumulate transformations */ L620: if (low > upp) { goto L700; } /* if */ if (notlas) { goto L660; } /* if */ for (i = low; i <= upp; ++i) { p = x * VECS(i, k ) + y * VECS(i, k + 1); VECS(i, k + 1) -= p * q; VECS(i, k) -= p; } /* for */ goto L700; L660: for (i = low; i <= upp; ++i) { p = x * VECS(i, k) + y * VECS(i, k + 1) + z * VECS(i, k + 2); VECS(i, k + 2) -= p * r; VECS(i, k + 1 ) -= p * q; VECS(i, k) -= p; } /* for */ L700: ;} /* for */ L720: goto L160; /* One root found */ L740: EIGVALS(en).re = x + t; H1(en, en) = EIGVALS(en).re; EIGVALS(en).im = 0.0; CNT(en) = its; en = na; goto L140; /* Two roots found */ L760: p = (y - x) / 2.0; q = p * p + w; z = sqrt(FABS(q)); x += t; H1(en, en) = x; H1(na, na) = y + t; CNT(en) = -its; CNT(na) = its; if (q < 0.0) { goto L840; } /* if */ /* Real pair */ if (p < 0.0) { z = p - z; } /* if */ if (p > 0.0) { z = p + z; } /* if */ EIGVALS(na).re = x + z; EIGVALS(en).re = EIGVALS(na).re; if (z != 0.0) { EIGVALS(en).re = x - w / z; } /* if */ EIGVALS(na).im = 0.0; EIGVALS(en).im = 0.0; x = H1(en, na); r = a02dbc(a02bac(x, z)); p = x / r; q = z / r; /* Row modification */ for (j = na; j <= n; ++j) { z = H1(na, j); H1(na, j) = q * z + p * H1(en, j); H1(en, j) = q * H1(en, j) - p * z; } /* for */ /* Column modification */ for (i = 1; i <= en; ++i) { z = H1(i, na); H1(i, na) = q * z + p * H1(i, en); H1(i, en) = q * H1(i, en) - p * z; } /* for */ /* Accumulate transformations */ for (i = low; i <= upp; ++i) { z = VECS(i ,na); VECS(i ,na) = q * z + p * VECS(i ,en); VECS(i ,en) = q * VECS(i ,en) - p * z; } /* for */ goto L860; /* Complex pair */ L840: EIGVALS(na).re = x + p; EIGVALS(en).re = x + p; EIGVALS(na).im = z; EIGVALS(en).im = -z; L860: en += -2; goto L140; /* All roots found now backsubstitute */ L880: if (norm == 0.0) { goto L1480; } /* if */ norm *= machep; /* Backsubstitution */ for (kk = 1; kk <= n; ++kk) { en = n + 1 - kk; p = EIGVALS(en).re; q = EIGVALS(en).im; na = en - 1; if (q != 0.0) { goto L1120; } /* if */ /* Real vector */ H1(en, en) = 1.0; if (na < 1) { goto L1340; } /* if */ for (ii = 1; ii <= na; ++ii) { i = na + 1 - ii; i1 = i - 1; w = H1(i, i) - p; r = H1(i, en); if (EIGVALS(i).im >= 0.0) { goto L900; } /* if */ z = w; s = r; goto L1100; L900: if (EIGVALS(i).im > 0.0) { goto L1020; } /* if */ /* Modification to stop overflow */ if (w != 0.0) { goto L940; } /* if */ if ((FABS(r)) < norm * 10.0) { goto L960; } /* if */ r = -r; for (j = 1; j <= en; ++j) { H1(j, en) *= norm; } /* for */ goto L980; L940: r = -r / w; goto L980; L960: r = -r / norm; L980: H1(i, en) = r; if (i1 == 0) { goto L1100; } /* if */ for (j = 1; j <= i1; ++j) { H1(j, en) += H1(j, i) * r; } /* for */ goto L1100; /* Solve real equations */ L1020: x = H1(i, i + 1 ); y = H1(i + 1, i ); q = (EIGVALS(i).re - p) * (EIGVALS(i).re - p) + EIGVALS(i).im * EIGVALS(i).im; t = (x * s - z * r) / q; H1(i, en) = t; if ((FABS(x)) > (FABS(z))) { goto L1040; } /* if */ r = (-s - y * t) / z; goto L1060; L1040: r = (-r - w * t) / x; L1060: H1(i+1, en) = r; if (i1 == 0) { goto L1100; } /* if */ for (j = 1; j <= i1; ++j) { H1(j, en) = H1(j, en) + H1(j, i + 1 ) * r + H1(j, i) * t; } /* for */ L1100: ;} /* for */ /* End real vector */ goto L1340; L1120: if (q > 0.0) { goto L1340; } /* if */ /* Complex vector associated with lambda=p-i*q */ if (FABS(H1(en, na)) <= FABS(H1(na, en))) { goto L1140; } /* if */ r = q / H1(en, na); s = -(H1(en, en) - p) / H1(en, na); goto L1160; L1140: tmpc1 = a02bac(0.0, -H1(na, en)); tmpc2 = a02bac(H1(na, na) - p, q); tmpc3 = a02cdc(tmpc1, tmpc2); r=a02bbc(tmpc3); s=a02bcc(tmpc3); L1160: H1(en, na) = 0.0; H1(en, en) = 1.0; H1(na, na) = r; H1(na, en) = s; if (na < 2) { goto L1340; } /* if */ na1 = na - 1; for (j = 1; j <= na1; ++j) { H1(j, en) += H1(j, na) * s; H1(j, na) *= r; } /* for */ for (ii = 1; ii <= na1; ++ii) { i = na1 + 1 - ii; i1 = i - 1; w = H1(i, i) - p; ra = H1(i, na); sa = H1(i, en); if (EIGVALS(i).im >= 0.0) { goto L1200; } /* if */ z = w; r = ra; s = sa; goto L1320; L1200: if (EIGVALS(i).im == 0.0) { goto L1280; } /* if */ /* Solve complex equations */ x = H1(i, i + 1); y = H1(i+1,i); vr = (EIGVALS(i).re - p) * (EIGVALS(i).re - p) + EIGVALS(i).im * EIGVALS(i).im - q * q; vi = (EIGVALS(i).re - p) * 2.0 * q; if (vr == 0.0 && vi == 0.0) { vr = machep * norm * ((FABS(w)) + (FABS(q)) + (FABS(x)) + (FABS(y)) + ( FABS(z))); } /* if */ tmpc1 = a02bac(x*r-z*ra+q*sa, x*s-z*sa-q*ra); tmpc2 = a02bac(vr,vi ); tmpc3 = a02cdc(tmpc1, tmpc2); t=a02bbc(tmpc3); u=a02bcc(tmpc3); if ((FABS(x)) <= (FABS(z)) + (FABS(q))) { goto L1220; } /* if */ r = (-ra - w * t + q * u) / x; s = (-sa - w * u - q * t) / x; goto L1240; L1220: tmpc1 = a02bac(-r-y*t, -s-y*u); tmpc2 = a02bac(z,q ); tmpc3 = a02cdc(tmpc1, tmpc2); r=a02bbc(tmpc3); s=a02bcc(tmpc3); L1240: H1(i, na) = t; H1(i, en) = u; H1(i+1, na) = r; H1(i+1, en) = s; if (i1 == 0) { goto L1320; } /* if */ for (j = 1; j <= i1; ++j) { H1(j, na) = H1(j, na) + H1(j, i + 1)* r + H1(j, i) * t; H1(j, en) = H1(j, en) + H1(j, i + 1) * s + H1(j, i) * u; } /* for */ goto L1320; L1280: tmpc1 = a02bac(-ra,-sa); tmpc2 = a02bac(w,q ); tmpc3 = a02cdc(tmpc1, tmpc2); r=a02bbc(tmpc3); s=a02bcc(tmpc3); H1(i, na) = r; H1(i, en) = s; if (i1 == 0) { goto L1320; } /* if */ for (j = 1; j <= i1; ++j) { H1(j, na) += H1(j, i) * r; H1(j, en) += H1(j, i) * s; } /* for */ L1320: ;} /* for */ /* End complex vector */ L1340: ;} /* for */ /* End backsubstitution */ /* Vectors of isolated roots */ low1 = low - 1; upp1 = upp + 1; for (j = 1; j <= n; ++j) { m = MIN(low1, j); if (m < 1) { goto L1380; } /* if */ for (i = 1; i <= m; ++i) { VECS(i ,j) = H1(i, j); } /* for */ L1380: if (upp1 > j) { goto L1420; } /* if */ for (i = upp1; i <= j; ++i) { VECS(i ,j) = H1(i, j); } /* for */ L1420: ;} /* for */ /* * Multiply by transformation matrix to give * vectors of original full matrix */ for (jj = low; jj <= n; ++jj) { j = low + n - jj; m = MIN(upp, j); for (i = low; i <= upp; ++i) { VECS(i ,j) = VECS(i ,m) * H1(m, j); } /* for */ --m; if (m + 1 >= low) { f06pac(NoTranspose, upp-low+1, m-low+1, 1.0, &VECS(low ,low), tdvecs, &H1(low, j), tdh, 1.0, &VECS(low ,j), tdvecs); } /* if */ } /* for */ L1480: /*if (fail && fail->code==NE_NOERROR) printf("In f02aqc: NE_NOERROR\n");*/ /*p01acc((char *)0, NE_NOERROR, "f02aqc", fail);*/ return ; L1500: Vsprintf(buf,"NE_TOO_MANY_ITERATIONS:\n\ More than %1ld iterations are required to isolate all the eigenvalues" ,30*n); printf("In %s: %s\n","f02aqc",buf); /*p01acc(buf, NE_TOO_MANY_ITERATIONS,"f02aqc",fail);*/ return ; } /* f02aqc_ */ /*------------------------------------------------------------------*/ /* New Routine */ /*------------------------------------------------------------------*/ /* f02agct.c * * Copyright 1994 Numerical Algorithms Group. * * NAG C Library * * Derived from: NAG Fortran Library F02AGF (Last Revision Mark 14A) * * Eigenvalues and eigenvectors of real unsymmetric matrix. * * C version: Shah Datardina NAG Ltd., Oxford, U.K. * * Mark 4 re-issue, 1994. * */ #ifdef NAG_PROTO void f02agc(long n, double a[], long tda, Complex r[], Complex v[], long tdv, long iter[], NagError *fail) #else void f02agc(n, a, tda, r, v,tdv, iter, fail) double a[]; long tda, n; Complex r[], v[]; long tdv; long iter[]; NagError *fail; #endif { /* Local variables */ double d__1, d__2; double term, c, d; double *vreal=(double *)0, *rreal=(double *)0; double machep, max_, sum; long c__1 = 1; long i, j, k, l; long ib; char buf[NAG_ERROR_BUF_LEN]; NagError local_fail; /* Function Body */ #define ITER(I) iter[(I)-1] #define RI(I) r[(I)-1].im #define RR(I) r[(I)-1].re #define RREAL(I) rreal[(I)-1].re #define VI(I,J) v[((I)-1) * ( tdv) + ((J)-1)].im #define VR(I,J) v[((I)-1) * ( tdv) + ((J)-1)].re #define VREAL(I,J) vreal[((I)-1) * ( n) + ((J)-1)] #define A4(I,J) a[((I)-1) * ( tda) + ((J)-1)] machep = X02AJC; ib = X02BHC; INIT_FAIL(local_fail); if (fail) fail->code = NE_NOERROR; if (n < 1) { Vsprintf(buf, "NE_INT_ARG_LT:\n\ On entry, %s must not be less than %s: %s = %ld.", "n", "1", "n", n); printf("In %s: %s\n","f02agc",buf); /*p01acc(buf, NE_INT_ARG_LT, "f02agc", fail);*/ return; } if (tda < n) { Vsprintf(buf, "NE_2_INT_ARG_LT:\n\ On entry %s = %2ld while %s = %2ld.\n\ These parameters must satisfy %s >= %s", "tda", tda, "n", n, "tda", "n"); printf("In %s: %s\n","f02agc",buf); /* p01acc(buf, NE_2_INT_ARG_LT, "f02agc", fail);*/ return; } if (tdv < n) { Vsprintf(buf,"NE_2_INT_ARG_LT:\n\ On entry %s = %2ld while %s = %2ld.\n\ These parameters must satisfy %s >= %s", "tdv", tdv, "n", n, "tdv", "n"); printf("In %s: %s\n","f02agc",buf); /* p01acc(buf, NE_2_INT_ARG_LT, "f02agc", fail);*/ return; } if (!(vreal= NAG_ALLOC(n*n, double))) { printf("In %s: %s\n","f02agc","NE_ALLOC_FAIL:\n\ Memory allocation failed."); /* p01acc((char *)0, NE_ALLOC_FAIL, "f02agc", fail);*/ return; } if (!(rreal= NAG_ALLOC(n, double))) { printf("In %s: %s\n","f02agc","NE_ALLOC_FAIL:\n\ Memory allocation failed."); /*p01acc((char *)0, NE_ALLOC_FAIL, "f02agc", fail);*/ return; } /* write_mat_r("In f02agf input a", 'c', a, tda, n*n, 0,1,n); */ f01atc(n, ib, a, tda, &k, &l, rreal); /* write_mat_r("On output from f01atf, a", 'c', a, tda, n*n, 0,1,n); */ /* write_mat_r("On output from f01atf, rreal", 'f',rreal, n, n, 0,1,1); */ f01akc(n, k, l, a, tda, iter); /* write_mat_r("On output from f01akf, a", 'c', a, tda, n*n, 0,1,n); */ /* write_mat_i("On output from f01akf, iter", 'f',iter, n, n, 0,1,1); */ f01apc(n, k, l, iter, a, tda, vreal, n); /* write_mat_r("On output from f01apf, vreal", 'c', vreal, n, n*n, 0,1,n); */ f01auc(n, k, l, n, rreal, vreal, n); /* write_mat_r("On output from f01auf, vreal", 'c', vreal, n, n*n, 0,1,n); */ f02aqc(n, c__1, n, machep, a, tda, vreal, n, r, iter, &local_fail); for(i=1; i<=n; i++) { rreal[i-1] = RR(i); } for(i=1; i<=n; i++) for(j=1; j<=n; j++) { VR(i,j)=VREAL(i,j); } /* write_mat_r("On output from f02aqf, vreal", 'c', vreal, n, n*n, 0,1,n); */ /* write_mat_r("On output from f02aqf, rreal", 'f',rreal, n, n, 0,1,1); */ /* write_mat_i("On output from f02aqf, iter", 'f',iter, n, n, 0,1,1); */ if (local_fail.code != NE_NOERROR) { Vsprintf(buf, "NE_TOO_MANY_ITERATIONS:\n\ More than %1ld iterations are required to isolate all the eigenvalues", 30*n); printf("In %s: %s\n","f02aqc",buf); /*p01acc(buf, NE_TOO_MANY_ITERATIONS,"f02agc",fail);*/ return; } for (i = 1; i <= n; ++i) { if (RI(i) == 0.0) { goto L60; } if (RI(i) > 0.0) { goto L100; } for (j = 1; j <= n; ++j) { VR(j,i) = VR(j,i-1); VI(j,i) = -VI(j,i-1); /* L40: */ } goto L140; L60: for (j = 1; j <= n; ++j) { VI(j,i) = 0.0; /* L80: */ } goto L140; L100: for (j = 1; j <= n; ++j) { VI(j,i) = VR(j,i+1); /* L120: */ } L140: ; } for (i = 1; i <= n; ++i) { sum = 0.0; max_ = 0.0; for (j = 1; j <= n; ++j) { if ((d__1 = VR(j,i), ABS(d__1)) <= max_) { goto L160; } max_ = (d__1 = VR(j,i), ABS(d__1)); L160: if ((d__1 = VI(j,i), ABS(d__1)) <= max_) { goto L180; } max_ = (d__1 = VI(j,i), ABS(d__1)); L180: ; } for (j = 1; j <= n; ++j) { VR(j,i) /= max_; VI(j,i) /= max_; /* L200: */ } max_ = 0.0; for (j = 1; j <= n; ++j) { d__1 = VR(j,i); d__2 = VI(j,i); term = d__1 * d__1 + d__2 * d__2; sum += term; if (term <= max_) { goto L220; } max_ = term; c = VR(j,i); d = -VI(j,i); L220: /* L240: */ ; } d__1 = c; d__2 = d; sum *= d__1 * d__1 + d__2 * d__2; sum = sqrt(sum); for (j = 1; j <= n; ++j) { term = VR(j,i); VR(j,i) = (VR(j,i) * c - VI(j,i) * d) / sum; VI(j,i) = (d * term + c * VI(j,i)) / sum; /* L260: */ } /* L280: */ } //NAG_FREE(rreal); //NAG_FREE(vreal); /*if (fail && fail->code==NE_NOERROR) p01acc((char *)0, NE_NOERROR, "f02agc", fail);*/ return; } /* f02agc */