Numerical Recipes Forum  

Go Back   Numerical Recipes Forum > Obsolete Editions Forum > Methods: Chapters 2, 11, and 18

Reply
 
Thread Tools Display Modes
  #1  
Old 04-02-2012, 11:35 PM
chatterjee chatterjee is offline
Registered User
 
Join Date: Mar 2012
Posts: 2
Problem using LUDCMP

Dear all,

I have been using NRCP recently.I am not familiar with the NR functions.As far as I could make out the functions from NRCP book,I tried that much.I'm trying to use LUDCMP function perform an LU Decomposition of a matrix.

Prototype Declaration

void ludcmp(float **a,int n,int *indx,float *d); (line 5)

Within Main

float **a,*d;
int *indx;
int i,j;

a = matrix(1,4,1,4);
indx = vector(1,4); (line 25)
d = vector(1,4);

a[4][4] = {{5.0,7.0, 4.0,11.0},{1.0,6.0,-19.0,1.0},{17.0,-3.0,9.1,21.0},{19.0,-6.0,10.5,32.0}}; (line 28)
d[4] = {3.40,-7.0,62.0,9.0}; (line 29)


ludcmp(a,4,indx,&d); (line 34)

...and then print the matrix accordingly.

But I get the following errors :

lu.c: In function ?main?:
lu.c:25:6: warning: assignment from incompatible pointer type [enabled by default]
lu.c:28:11: error: expected expression before ?{? token
lu.c:29:8: error: expected expression before ?{? token
lu.c:34:1: warning: passing argument 4 of ?ludcmp? from incompatible pointer type [enabled by default]
lu.c:5:6: note: expected ?float *? but argument is of type ?float **?

Please help!
Reply With Quote
  #2  
Old 04-03-2012, 10:28 AM
davekw7x davekw7x is offline
Registered User
 
Join Date: Jan 2008
Posts: 453
Quote:
Originally Posted by chatterjee View Post
...NRCP book...
I don't know anything about the NRCP book. The Numerical Recipes Version 2 C book is available for download from the nr.com web site. I suggest you read that one.

Anyhow...

The very old Version 2 of Numerical Recipes C library functions were more-or-less transmogrified from the older Fortran programs, and so the NR Version 2 distribution used vectors and matrices that were 1-based, not 0-based. That is, the NR functions expect that an NR vector with size N will have index values 1, 2, ... , N-1, and an NR matrix with size NxN will have index values 1 through N-1. (See Footnote.)

Furthermore...

There is a function, convert_matrix() (in nrutil.c) that creates a matrix and initializes it from a C array (1-D or 2-D) that contains matrix values. On the other hand, I usually just declared and initialized a 2-D array in C and then created a NR matrix and copied array elements into the matrix manually. Keeps me grounded.

I will note here that, in C (or C++), you can not assign multiple values to multiple elements of an array (of any dimension) by using the array name in an assignment statement. I mean, you can initialize arrays as I show in the following example, but the comma-separated-list-enclosed-in-braces notation only works when the array is declared; that notation can not (that's not) be used in any kind of assignment statement. Period.

An example that illustrates these points:

Code:
/*
    NR Version 2 C example showing invocation of ludcmp()

    davekw7x
*/
#include "nr.h"
#include "nrutil.h"

/* Size of matrix is NP x NP */
const int NP = 4;

int main()
{
    float **a; /* Input and Output: The matrix                              */
    int *indx; /* Output: Keeps track of swapped rows                       */
    float d;   /* Output: Indicates whether number of swaps was even or odd */

    int i, j;  /* Loop index variables                                      */

    /* C-style array with initialization */
    float a_array[4][4] = {
        { 5.0,  7.0,   4.0, 11.0},
        { 1.0,  6.0, -19.0,  1.0},
        {17.0, -3.0,   9.1, 21.0},
        {19.0, -6.0,  10.5, 32.0}
    };

    /* Create NR objects */
    a = matrix(1, NP, 1, NP);
    indx = ivector(1, NP);

    /* Copy C-style array values into the NR object */
    for (i = 1; i <= NP; i++) {
        for (j = 1; j <= NP; j++) {
            a[i][j] = a_array[i-1][j-1];
        }
    }

    /* Print the original matrix */
    printf("Initially:  [A]\n");
    for (i = 1; i <= NP; i++) {
        for (j = 1; j <= NP; j++) {
            printf("  %+e", a[i][j]);
        }
        printf("\n");
    }
    printf("\n");

    /* Do the decomposition */
    ludcmp(a, NP, indx, &d);

    /* Print the result of the decomposition */
    printf("After ludcmp: [A]\n");
    for (i = 1; i <= NP; i++) {
        for (j = 1; j <= NP; j++) {
            printf("  %+e", a[i][j]);
        }
        printf("\n");
    }
    printf("\n");

    /*
       It's considered a good habit to get into:
       Deallocate storage for all dynamic objects
    */
    free_ivector(indx, 1, NP);
    free_matrix(a, 1, NP, 1, NP);

    return 0;
}
(Note that the vector() and matrix() functions create floating point objects. If you want a vector of ints, use ivector(),as I show above.)



Output using Numerical Recipes Version 2.11 official distribution and after compiling with GNU gcc version 4.1.2 on my Linux workstation:
Code:
Initially:  [A]
  +5.000000e+00  +7.000000e+00  +4.000000e+00  +1.100000e+01
  +1.000000e+00  +6.000000e+00  -1.900000e+01  +1.000000e+00
  +1.700000e+01  -3.000000e+00  +9.100000e+00  +2.100000e+01
  +1.900000e+01  -6.000000e+00  +1.050000e+01  +3.200000e+01

After ludcmp: [A]
  +1.700000e+01  -3.000000e+00  +9.100000e+00  +2.100000e+01
  +2.941177e-01  +7.882353e+00  +1.323529e+00  +4.823529e+00
  +5.882353e-02  +7.835821e-01  -2.057239e+01  -4.014925e+00
  +1.117647e+00  -3.358209e-01  -3.761742e-02  +9.998222e+00
Regards,

Dave

Footnote:
Don't like the 1-based array stuff? Well, neither does anyone else (or at least any other C/C++ programmers that I know). That's why we welcomed the release of NR Version 3 a few years ago. The NR Version 2 C++ distribution used 0-based arrays also, but C programmers hated to have to figure out some of the things that you had to do to use NR2 C++, and C++ programmers hated the fact that the NR2 C++ stuff didn't take advantage of very much of the power of C++.

In my opinion, NR3 is much more acceptable (mainly since the hardcore C guys seem to have accepted promotions to management positions rather than actually learn some significant C++). I mean, in addition to the real advantage of having a fair number of new algorithms and improved implementations of some of the old ones, NR3 programs use templates and functors and stuff like that.

Of course, unlike the apparently free NRCP stuff, NR3 is not free. Oh, well...

Last edited by davekw7x; 04-04-2012 at 12:01 AM.
Reply With Quote
  #3  
Old 04-07-2012, 01:13 AM
chatterjee chatterjee is offline
Registered User
 
Join Date: Mar 2012
Posts: 2
Thanks Dave for your detailed reply.Being new to NR,I messed up with the array index.Now with the help of your program I have solved the matrix with LUBKSB and found out the inverse as well. The "NRCP book" I was talking about is actually the "Numerical Recipes in C" book.

Away from the discussion,I seek your opinion on the other free libraries available for C like GNU Scientific Library.Have you used it? How does it differ from NR3 as far as the number of functions and usage is concerned?
Reply With Quote
Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off

Forum Jump


All times are GMT -5. The time now is 08:33 PM.


Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.