Inverse Matrix Implementation Problem...

Inverse Matrix Implementation Problem...

[W L A Au]

[-- Attachment #1: Type: text/plain, Size: 2151 bytes --]

Dear Sir/Madam,

Please help me. I've encounter some problems when I try to implement the
inverse matrix function by using the Gauss-Jordan Elimination routine,
which is provided by the Numerical Recipes in C, Chapter 2.1

http://www.ulib.org/webRoot/Books/Numerical_Recipes/bookcpdf.html

What exactly I want to do is:

First, put in the values of a 2D matrix, A

                3 5 0
for example, A= 6 1 4
                1 1 2

Then, add an equal dimension 2D identity matrix to A

   3 5 0 1 0 0
=> 6 1 4 0 1 0
   1 1 2 0 0 1
   
After that, do row elementary operation (eg. addition, multiplication, swap
row...) in order to form an identity matrix on the left

   1 0 0 1/23  5/23  -10/23
=> 0 1 0 4/23  -3/23 6/23
   0 0 1 -5/46 -1/23 27/46

When the left hand side become an identity matrix, the right hand side will
become the inverse of A. So, I use the gaussj() function
provided by the book to achieve this. However, when I use the gcc compiler
command:

gcc test2.c -o test2

I've got the following error message:

Undefined			first referenced
 symbol  			    in file
free_ivector                        /var/tmp/cc4v8fRq.o
ivector                             /var/tmp/cc4v8fRq.o
nrerror                             /var/tmp/cc4v8fRq.o
ld: fatal: Symbol referencing errors. No output written to test2
collect2: ld returned 1 exit status

It seems I can not include the nrutil.h header file. But I've put it in the
same subdirectory with test2.c already. How this can be fix?

Also, the 2D matrix, a. Is already declared as "float **a" in the routine.
I've try to use it by substitude values in, like

a[0][0] = 3;   a[1][0] = 5;   a[2][0] = 0;
a[0][1] = 6;   a[1][1] = 1;   a[2][1] = 4;
a[0][2] = 1;   a[1][2] = 1;   a[2][2] = 2;

in order to find the inverse matrix of a. But error massage come out as
well. I guess this is not right at all. How can I make it right?

Finally, I'm not sure did I use the gaussj() function correctly. Did I miss
something in my implementation?

I would like to thank you in advance for helping me. Thank a lot!

Yours Sincerely,
Alderick

PS. The test2.c and nrutil.h files are attached with this email.

[-- Attachment #2: test2.c --]
[-- Type: , Size: 1739 bytes --]

#include <stdlib.h>
#include <stdio.h>

#include <math.h>
#include "nrutil.h"

#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}

void gaussj(float **a, int n, float **b, int m)
{
    int *indxc,*indxr,*ipiv;
    int i,icol,irow,j,k,l,ll;
    float big,dum,pivinv,temp;

    indxc=ivector(1,n);
	indxr=ivector(1,n);
    ipiv=ivector(1,n);
    for (j=1;j<=n;j++) ipiv[j]=0;
    for (i=1;i<=n;i++) {
		big=0.0;
        for (j=1;j<=n;j++)
			if (ipiv[j] != 1)
				for (k=1;k<=n;k++) {
					if (ipiv[k] == 0) {
						if (fabs(a[j][k]) >= big) {
							big=fabs(a[j][k]);
							irow=j;
							icol=k;
						}
					} else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1");
				}
		++(ipiv[icol]);

        if (irow != icol) {
			for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
			for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
		}
		indxr[i]=irow;
		indxc[i]=icol;
		if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2");
		pivinv=1.0/a[icol][icol];
		a[icol][icol]=1.0;
		for (l=1;l<=n;l++) a[icol][l] *= pivinv;
		for (l=1;l<=m;l++) b[icol][l] *= pivinv;
		for (ll=1;ll<=n;ll++)
			if (ll != icol) {
				dum=a[ll][icol];
				a[ll][icol]=0.0;
				for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
				for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
			}
	}

	for (l=n;l>=1;l--) {
		if (indxr[l] != indxc[l])
			for (k=1;k<=n;k++)
				SWAP(a[k][indxr[l]],a[k][indxc[l]]);
	}
	free_ivector(ipiv,1,n);
	free_ivector(indxr,1,n);
	free_ivector(indxc,1,n);
}

main()
{
    float **a, **b;
    int n, m;

    a[0][0] = 3;   a[1][0] = 5;   a[2][0] = 0;
    a[0][1] = 6;   a[1][1] = 1;   a[2][1] = 4;
    a[0][2] = 1;   a[1][2] = 1;   a[2][2] = 2;

    gaussj(a, n, b, m);
}

[-- Attachment #3: nrutil.h --]
[-- Type: , Size: 2995 bytes --]

// Here is a listing of nrutil.h:

#ifndef _NR_UTILS_H_
#define _NR_UTILS_H_

static float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)


static double dsqrarg;
#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)


static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ? (dmaxarg1) : (dmaxarg2))


static double dminarg1,dminarg2;
#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ? (dminarg1) : (dminarg2))


static float maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))


static float minarg1,minarg2;
#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ? (minarg1) : (minarg2))


static long lmaxarg1,lmaxarg2;
#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ? (lmaxarg1) : (lmaxarg2))


static long lminarg1,lminarg2;
#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ? (lminarg1) : (lminarg2))


static int imaxarg1,imaxarg2;
#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ? (imaxarg1) : (imaxarg2))


static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ? (iminarg1) : (iminarg2))

#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

#if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */

void nrerror(char error_text[]);
float *vector(long nl, long nh);
int *ivector(long nl, long nh);
unsigned char *cvector(long nl, long nh);
unsigned long *lvector(long nl, long nh);
double *dvector(long nl, long nh);
float **matrix(long nrl, long nrh, long ncl, long nch);
double **dmatrix(long nrl, long nrh, long ncl, long nch);
int **imatrix(long nrl, long nrh, long ncl, long nch);
float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch, long newrl, long newcl);
float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_vector(float *v, long nl, long nh);
void free_ivector(int *v, long nl, long nh);
void free_cvector(unsigned char *v, long nl, long nh);
void free_lvector(unsigned long *v, long nl, long nh);
void free_dvector(double *v, long nl, long nh);
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh);

#else /* ANSI */
/* traditional - K&R */

void nrerror();
float *vector();
// Rest of traditional declarations are here on the diskette.

#endif /* ANSI */

#endif /* _NR_UTILS_H_ */

Re: Inverse Matrix Implementation Problem...

[zhenggen]

W L A Au wrote:

>Dear Sir/Madam,
>
>Please help me. I've encounter some problems when I try to implement the
>inverse matrix function by using the Gauss-Jordan Elimination routine,
>which is provided by the Numerical Recipes in C, Chapter 2.1
>
>http://www.ulib.org/webRoot/Books/Numerical_Recipes/bookcpdf.html
>
>What exactly I want to do is:
>
There are two basic concepts which you should grasp before you can 
continue your work.
(1) The declaration of array(one-dimensional or multi-dimensional):
  "  float **a ;" will not declare an array, rather it declares a 
pointer. "a" is a pointer which can point to an array of float pointers.
 The simplest way to declare a 3 by 3 array is :    float a[3][3];
(2) header files and library files.
   In a .C file, a "#include ..." statement only tells the compiler to 
find ,in this header file, the definition of symbols appeared in this C 
file . The compile stage will pass smoothly if you don't provide the 
respective object file, which is normally put   in an archived library 
file together with other object files. The object file, which is created 
from the relevant C file, is needed in link phase. If the archived 
library file is named libmath.a, then you have to tell ld to attach this 
library in the following way:
    gcc test2.c -o test2 -lmath

Please notice the notation. It should be "-lmath" ,not "-llibmath.a"
In your codes, there is a line:

#include "nrutil.h"

but no library is given explicitly in the command line. Surely you will 
fail.
  There are one exception. When you put line like "#include <stdlib.h>" 
in your c file, the relavant object file, or library will be 
automatically attached to your execute file. It is called C runtime library.

    I hope this will help you.

Re: Inverse Matrix Implementation Problem...

[Edmund Green]

re:
> Please help me. I've encounter some problems when I try to implement the
> inverse matrix function by using the Gauss-Jordan Elimination routine,
> which is provided by the Numerical Recipes in C, Chapter 2.1
...
 > However, when I use the gcc compiler command:
> gcc test2.c -o test2
> 
> I've got the following error message:
> 
> Undefined			first referenced
>  symbol  			    in file
> free_ivector                        /var/tmp/cc4v8fRq.o
> ivector                             /var/tmp/cc4v8fRq.o
> nrerror                             /var/tmp/cc4v8fRq.o
> ld: fatal: Symbol referencing errors. No output written to test2
> collect2: ld returned 1 exit status

These missing library routines are found in the nrutil.c file (see 
Appendix B of the book), to save having to duplicate them for each 
example in the book.

This makes your compilation more complex, you should really be using 
makefiles for anything with more than 1 source file (see "info make").
They are a very powerful but can also get very complicated.

However in this simple case, if you don't mind 'make' doing a few things 
behind your back with its built in implicit rules, create a file called 
"makefile" with the following single line in it

test2 : test2.o nrutil.o

and put in in a directory that also contains "test2.c", "nrutil.h" and 
"nrutil.c", then run the command "make".


Edmund.