A Scientific Computation Library

Ky Van Ha

1. Introduction

Scientific Computation becomes more and more important to our students. A Scientific Computation Library (SCL) that they can use in their applications is clearly of interest. I started a project six months ago. The goal is to develop a scientific computation library with C++ which handles some problems in Applied Mathematics, Numerical Analysis, Graphics, Signal Processing, Image Processing, Wavelet Analysis, Neural Network, etc.

The Library will be updated and expanded each year, and a document will be made to report on the new results. In order that such a library can be succesfully developed,many people, students and colleagues, must be hopefully involved. The Library that is on the net now is really an extension of the Graphics Library I developed last year, see [1], it contains more than 150 functions. The Library has two versions, one for DOS and one for Windows.

The Library has been contructed and tested with Borland C++4.5, using the standard API functions. Thus, with a little change, it can be used for other compilators as Microsoft C++ as well.

In this paper, I would like to give some examples to show how one can use that library in Windows Applications, see [1] and [2] for more detail.

2. Compile an application with SCL

There are fives prototypes files one can use. These files are on the drive "I" under the directory "bc45\include\kvh" :

kvh_def.h   :  Global definitions and macros 
kvh_wmat.h  :  Mathematical routines 
kvh_wgrp.h  :  Graphics routines 
kvh_wnum.h  :  Numerical routines 
kvh_wsgn.h  :  Signal and wavelet routines

In order to compile your application with the SCL, you must :

include one or many prototype files in your source code file, for example :
```
#include <kvh\kvh_wsgn.h> 
#include <kvh\kvh_wgrp.h>
```
include the following "lib" file in your "*.ide" file :
i:\bc45\lib\kvh\kvh_wlib[.lib]

When you write a scientific program, a book like " Numerical Recipe in C" (see [3]) is of course very helpful. But the book is not easy for everybody. With the SCL, you can just call the functions that are defined in the prototype files from your programs.

In the following sections, some examples about graphics, mathematics and signal processing will be given using the scientific computation library. These examples demostrate also how one can write a Windows program with Borland C++4.5

3. Using the numerical routines

The following program will solve the linear system Mx = b using LU factorization.


#include <owl\applicat.h> 
#include <owl\framewin.h> 
#include <owl\dc.h> 
#include <kvh\kvh_wnum.h> 
 
// Function to solve a linear system 
// using LU factorization 

void Linear_System(float *x) 
{ 
    int n = 4, i, j; 
    float **A, **u, **l; 
    float M[][4] = {{ 3, -13, 9, 3}, 
                    {-6, 4, 1, -18}, 
                    { 6, -2, 2, 4 }, 
                    {12, -8, 6, 10}}; 
    float b[] = {-19, -34, 16, 26}; 
 
    A = dmatrix(n, n); 
    u = dmatrix(n, n); 
    l = dmatrix(n, n);
 
    for (i = 0; i < n; i++) 
        for (j = 0; j < n; j++)
	    A[i][j] = M[i][j]; 
 
    LU_fac(A, u, l, n); 
    LU_solve(l, u, b, x, n); 
 
    delete A; 
    delete u; 
    delete l; 
} 
 
class TMyWindow : public TFrameWindow 
{ 
    public: 
        TMyWindow(TWindow *Parent, LPSTR Title); 
 
    protected: 
	void Paint(TDC& DC, BOOL, TRect&); 
	void EvSize(UINT sizeType, TSize& size); 
 
	DECLARE_RESPONSE_TABLE(TMyWindow); 
}; 
 
DEFINE_RESPONSE_TABLE1(TMyWindow, TFrameWindow) 
    EV_WM_SIZE, 
END_RESPONSE_TABLE; 

TMyWindow::TMyWindow( TWindow* Parent,
                      LPSTR Title ) : 
    TFrameWindow(Parent, Title) 
{ 
    Attr.W = Attr.H = 400; 
}; 

void TMyWindow::EvSize(UINT sizeType, TSize& size) 
{ 
    TFrameWindow::EvSize(sizeType, size); 
    Invalidate(); 
} 
 
void TMyWindow::Paint(TDC& DC, BOOL, TRect&) 
{ 
    float x[4]; 
    char Msg[100]; 
 
    Linear_System(x); // Solve Mx = b, 
                      // solutions are stored in x 
 
    sprintf(Msg,"%.3f  %.3f  %.3f  %.3f\n", 
            x[0], x[1], x[2], x[3]); 
 
    DC.SetTextColor(TColor::LtBlue); 
    DC.TextOut(10, 50, Msg, lstrlen(Msg)); 
}; 
 
// class TApp 
class TApp : public TApplication 
{ 
    public: 
        TApp(): TApplication() {} 
	void InitMainWindow(); 
}; 
 
// This function creates the 
// application's main window. 
void TApp::InitMainWindow() 
{ 
    TFrameWindow *wndw = new TMyWindow(0, "Test"); 
    SetMainWindow(wndw); 
} 
 
int OwlMain(int, char* [] ) 
{ 
    return TApp().Run(); 
}

The program is a standard Windows program. The Paint function will write the results in the windows after it called the function "Linear_System(x)" that we define in the beginning of the program. That function calls two numerical routines defined in the SCL:


LU_fac(A, u, l, n);
// Compute A = LU, LU factorization
LU_solve(l, u, b, x, n);
// Solve Ax = b by LUx = b

The solutions are stored in the array x. Write it to a string Msg, and use DC.TextOut function to write this string to Windows.

The result is 3.000 1.000 -2.000 1.000

4. Using the graphics routines

In the next example, we would like to interpolate a data set using B-spline technique. You can define a class, say TBspline, that does everything for you. It can generate a B-spline curve, interpolate a data set, approximate a function, compute derivative, integral by calling the very technical routines in the SCL.


 
#if !defined (_TBSPLINE_H) 
#define _TBSPLINE_H 
 
#include  
#include  
#include  
 
// TBspline Class 
 
class TBspline 
{ 
private: 
    int NumControlP;    // n 
    int NumPnt; 
    int Bspline_Order;  // k 
    float *t;		// Knot vector 
    float *b; 
    RPoint *ControlP; 
    RPoint *Rpt; 
    void Gen_RPt(float *x, float *y); 
public: 
    TBspline(int numpnt, int numcon, int order);
    ~TBspline(); 
    void Gen_Bspl_Curve(RPoint *CP);
    void Bspl_Inter(float *x, float *y);
    void Bspl_Appro(float (*f)(float),
                           float x0, float xn);
    float Bspl_Deriv(float (*f)(float),
                            float x0, float xn, 
			    float xx);
    float Bspl_Integ(float (*f)(float), 
                            float x0, float xn);
    void Display(HDC hDC, COLORREF);
    void Display_ControlNet(HDC hDC, COLORREF);
}; 
 
#endif

Let us look at an example. To implement the function Bspl_Inter(x, y), where x, y are two arrays which contain the coordinates of the data points, we will call the B-spline routines from the SCL.

void TBspline::Bspl_Inter(float *x, float *y) 
{ 
    // Compute the knot vector	 
    knot(NumControlP, Bspline_Order, t, x);

    // Interpolate data points 
    bspl_intp(x, y, t, NumControlP, 
              Bspline_Order);

    // Generate B-spline curve 
    Gen_RPt(x, y); 
}

In order to use the TBspline class, we have to write a function in our application, say DrawCurve, which will interpolate data points and draw a B-spline curve through these points. Here is an example:

 
void TMyWindow::DrawCurve(HDC hDC, RECT rect) 
{ 
    RPoint Left = { -1.5, -10, 0}; 
    RPoint Right = { 1.5, 15, 0}; 

    _RCoord(hDC, rect, Left, Right); 
    _Display_Axis(hDC, -1, -8, 1, 13, 0.5, 2, 
                  TRUE, RGB(0,0,0)); 

    _Set_Cross(hDC, n, x, y, RGB(255, 0, 0)); 
    // Set Cross at the data points

    Spl->Bspl_Inter(x, y); 
    Spl->Display(hDC, RGB(100, 10, 100)); 
};

The result is shown in the figure 1.

Figure 1: Interpolation using B-spline method

5. Using the signal processing routines

To design a filter is a common problem in signal processing. In the following esample, we will design a Lowpass Filter based on Kaiser Window. The program will compute the impulse response h(n) of the filter with the following design specifications:

_p = 0.2, A_s = 40dB and ₁ = ₂

(see [4], kap. 3). The program calls first the routine:

Kaiser_Win_Parm(40, 0.2*PI, 0.3*PI, &N, 
                &beta, &delta2);

that computes the Kaiser window parameters, and then use the parameter "beta" and the modified zeroth-order Bessel function to compute the Kaiser window function v(n). And finally it computes the impulse response h(n). The following codes do that jobb.

int N; 
float beta, delta2;
Kaiser_Win_Parm(40, 0.2*PI, 0.3*PI, &N, 
                &beta, &delta2);   
float B0 = Bess0(beta, 20); 
float *x = new float[N]; 
int N1 = N >> 1; 
  
float v, u; 
h[0] = 0.25; 
 
for (int i = 1; i <= N1; i++) 
{ 
    u = 1 - i*i/(N1*N1); 
    u = sqrt(u)*beta; 
    v = Bess0(u, 20)/B0; 
    h[i] = v*sin(PI*i/4)/(PI*i); 
}

The impulse response of the filter is displayed in figure 2.

Figure 2: Filter Design using Kaiser window

References

Ky Van Ha. Graphics Library for Windows, Hiø 1994
Ky Van Ha, Scientific Computation Library 1995, Hiø 1995
William H. Press et all, Numerical Recipe in C, Cambridge Uni. Press, 1992
P.P Vaidyanathan, Multirate systems and filter banks, Prentice Hall, inc. 1993

HØit Nr. 2-96