Page 1 of 2 12 LastLast
Results 1 to 10 of 17

Thread: C --> Timing a Matrix Inversion

  1. #1
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152

    C --> Timing a Matrix Inversion

    [font=courier new][size=8pt]Here is a C program which fills an N x N matrix with random doubles, and then,
    inverts it, using the Gauss-Jordan Method.

    It performs the inversion, and then calculates and writes an N-dimensional check
    vector, to an output file. If the check vector contains all '1's, then, the
    inversion is correct.

    I compiled it using Bloodshed Dev-C++, version 4.9.9.2.

    http://www.bloodshed.net/

    The compiled '.exe' file, is only 20 KB.

    On my computer, which was manufactured in 1945 (maybe a slight exaggeration), it
    takes in the range of, 33 to 37 seconds, to fill a 1000 x 1000 matrix with
    random doubles, and to invert it.

    I attached the source file below, in case anyone wants to try running the
    program. I called the attachment file, "main.txt". If you download it, then,
    change its name to "main.c". Put "main.c" into some folder. In the Bloodshed
    IDE, you have to make a new project (File --> New --> Project). You can call the
    project, anything you want. In the dialog box, choose "C", not, "C++"; and,
    "Console application". Then, choose the folder that contains "main.c", as the
    project folder. The IDE will then create and open a default file, called,
    "main.c". Close that file, without saving it. Then, add your version of "main.c"
    (from inside the project folder), to the project (Project --> Add to Project).
    Then, open, "main.c". Now, to run the program, just choose, "Execute --> Compile
    & Run" (F9). When you run it, a console window will open. The program will write
    its progress in the window, and will tell you when it finishes. After it
    finishes running, you can open the output file (by default, called,
    "matinv.txt"), in the IDE. That file gives the run's statistics (matrix size,
    elapsed time), and lists the calculated N-dimensional check vector. If the
    inversion is correct, then all of the N check vector values, should be very
    close to 1.


    Dan

    [code=C]//------------------------------------------------------------------------------------------------

    // A 'C' program to time matrix inversions.
    // Uses the Gauss-Jordan Method to invert an N x N matrix, filled with random doubles.
    // The results are written to a file; by default, "matinv.txt".

    //------------------------------------------------------------------------------------------------

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

    //------------------------------------------------------------------------------------------------

    // Set size of N x N matrix to be inverted.
    #define N 1000

    //------------------------------------------------------------------------------------------------

    const int first = 0;
    const int last = N - 1;

    // Name of outfile.
    const char outfilename[80] = "matinv.txt";

    //------------------------------------------------------------------------------------------------

    void invert(double original[N][N], double inversion[N][N]);
    void initializeall(double m1[N][N], double m2[N][N], double m3[N][N], double v1[N]);
    void setpivot(double m1[N][N], double m2[N][N], int col);
    void switchrows(double m1[N][N], int row1, int row2);
    void dividerow(double m1[N][N], int row, double divisor);
    void torowaddrow(double m1[N][N], int torow, int addrow, double addrowfactor);
    void matmatmult(double m1[N][N], double m2[N][N], double m3[N][N]);
    void matvecmult(double m1[N][N], double v1[N], double v2[N]);
    double getelapsedtime();
    void writeoutfile(double t, double v1[N]);

    //------------------------------------------------------------------------------------------------

    // Matrix to be inverted (filled with random doubles).
    double a0[N][N];

    // Copy of matrix to be inverted.
    double am[N][N];

    // Calculated inversion of a0.
    double ai[N][N];

    // Calculated identity matrix.
    double id[N][N];

    // va is multiplied by id, to give the check vector, vb.
    double va[N];

    // check vector.
    double vb[N];

    //------------------------------------------------------------------------------------------------

    int main(int argc, char *argv&#91;])
    {
    double t0, t1, tt;

    printf("Inverting matrix..\n");
    printf("\n");

    t0 = getelapsedtime();

    // a0 is filled with random doubles.
    // am is set equal to a0.
    // ai is set to the N x N identity matrix.
    // va is set to an N-dimensional vector, of all '1's.
    initializeall(a0, am, ai, va);

    // am becomes the N x N identity matrix.
    // ai becomes the (N x N) inversion of a0.
    invert(am, ai);

    t1 = getelapsedtime();
    tt = t1 - t0;

    printf("Finished inverting a %d by %d matrix.\n", N, N);
    printf("elapsed time = %6.3f seconds\n", tt);
    printf("\n");
    printf("Calculating check vector..\n");
    printf("\n");

    // Multiply the inverted matrix, ai, by the original matrix, a0, to give the calculated identity matrix, id.
    matmatmult(a0, ai, id);

    // Multiply the N-dimensional vector of all '1's, va, by the calculated N x N identity matrix, id, to give the check vector, vb.
    // If the inversion is correct, then, vb, should also be an N-dimensional vector, of all '1's.
    // vb, is written to the output file.
    matvecmult(id, va, vb);

    writeoutfile(tt, vb);

    printf("Finished calculating check vector.\n");
    printf("Output written to file, \"%s\".\n", outfilename);
    printf("\n");
    printf("Press 'Enter' to quit.\n");
    getchar();
    return 0;
    }

    //------------------------------------------------------------------------------------------------

    void invert(double original[N][N], double inversion[N][N])
    // Uses the Gauss-jordan Method.
    {
    int row, col;
    double divisor, factor;

    for(col = first; col <= last; col++)
    {
    setpivot(original, inversion, col);
    divisor = original[col][col];
    dividerow(original, col, divisor);
    dividerow(inversion, col, divisor);

    for(row = col + 1; row <= last; row++)
    {
    factor = -original[row][col];
    torowaddrow(original, row, col, factor);
    torowaddrow(inversion, row, col, factor);
    }
    }
    for(col = last; col >= 1; col--)
    for(row = col - 1; row >= first; row--)
    {
    factor = -original[row][col];
    torowaddrow(original, row, col, factor);
    torowaddrow(inversion, row, col, factor);
    }
    }

    //------------------------------------------------------------------------------------------------

    void initializeall(double m1[N][N], double m2[N][N], double m3[N][N], double v1[N])
    {
    int i, j;
    time_t t;


    srand(time(&t));

    for(i = first; i <= last; i++)
    for(j = first; j <= last; j++)
    {
    m1[i][j] = rand();
    m2[i][j] = m1[i][j];
    }

    for(i = first; i <= last; i++)
    {
    v1[i] = 1;
    for(j = first; j <= last; j++)
    {
    m3[i][j] = 0;
    if(i == j) m3[i][j] = 1;
    }
    }
    }

    //------------------------------------------------------------------------------------------------

    void setpivot(double m1[N][N], double m2[N][N], int col)
    {
    int row, maxrow;
    double max = 0;
    double test;

    for(row = col; row <= last; row++)
    {
    test = abs(m1[row][col]);
    if( test > max)
    {
    max = test;
    maxrow = row;
    }
    }
    if(maxrow != col)
    {
    switchrows(m1, col, maxrow);
    switchrows(m2, col, maxrow);
    }
    }

    //------------------------------------------------------------------------------------------------

    void switchrows(double m1[N][N], int row1, int row2)
    {
    double temp;
    int i;
    for(i = first; i <= last; i++)
    {
    temp = m1[row1][i];
    m1[row1][i] = m1[row2][i];
    m1[row2][i] = temp;
    }
    }

    //------------------------------------------------------------------------------------------------

    void dividerow(double m1[N][N], int row, double divisor)
    {
    int i;
    for(i = first; i <= last; i++) m1[row][i] /= divisor;
    }

    //------------------------------------------------------------------------------------------------

    void torowaddrow(double m1[N][N], int torow, int addrow, double addrowfactor)
    {
    int i;
    for(i = first; i <= last; i++) m1[torow][i] += m1[addrow][i] * addrowfactor;
    }

    //------------------------------------------------------------------------------------------------

    void matmatmult(double m1[N][N], double m2[N][N], double m3[N][N])
    {
    int i, j, k;
    for(i = 0; i <= last; i++)
    {
    for(j = 0; j <= last; j++)
    {
    m3[i][j] = 0;
    for(k = 0; k <= last; k++) m3[i][j] += m1[i][k] * m2[k][j];
    }
    }
    }

    //------------------------------------------------------------------------------------------------

    void matvecmult(double m1[N][N], double v1[N], double v2[N])
    {
    int i, j;
    for(i = 0; i <= last; i++)
    {
    v2[i] = 0;
    for(j = 0; j <= last; j++) v2[i] += m1[i][j] * v1[j];
    }
    }

    //------------------------------------------------------------------------------------------------

    double getelapsedtime()
    {
    time_t ticks;
    double et;
    ticks = clock();
    et = (double)ticks / CLOCKS_PER_SEC;
    return et;
    }

    //------------------------------------------------------------------------------------------------

    void writeoutfile(double t, double v1[N])
    {
    int i;
    FILE *outfile;

    outfile = fopen(outfilename, "w");
    fprintf(outfile, "Inversion of a %d by %d matrix,\n", N, N);
    fprintf(outfile, "filled with uniformly distributed random doubles,\n");
    fprintf(outfile, "in the range of, [0, 32767].\n");
    fprintf(outfile, "\n");
    fprintf(outfile, "elapsed time = %6.3f seconds\n", t);
    fprintf(outfile, "\n");
    fprintf(outfile, "If the inversion was calculated correctly, then every\n");
    fprintf(outfile, "value shown below should be very close to 1.\n");
    fprintf(outfile, "\n");

    for(i = first; i <= last; i++) fprintf(outfile, "%5i %19.16f\n", i + 1, v1[i]);

    fclose(outfile);
    }

    //------------------------------------------------------------------------------------------------
    [/code]
    Attached Files Attached Files
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  2. #2
    thinBasic MVPs kryton9's Avatar
    Join Date
    Nov 2006
    Location
    Naples, Florida & Duluth, Georgia
    Age
    67
    Posts
    3,869
    Rep Power
    404

    Re: C --> Timing a Matrix Inversion

    dev-c++ is very nice. If you want to make GUI c++ apps you might want to look at wxDev-c++
    It is dev-c++ with wxWidgets packaged in with nice wizards to get you started.
    http://wxdsgn.sourceforge.net/
    Acer Notebook: Win 10 Home 64 Bit, Core i7-4702MQ @ 2.2Ghz, 12 GB RAM, nVidia GTX 760M and Intel HD 4600
    Raspberry Pi 3: Raspbian OS use for Home Samba Server and Test HTTP Server

  3. #3

    Re: C --> Timing a Matrix Inversion

    running Windows XP(32-bit) under VMware Fusion on my mac pro it takes about 21 seconds, the same program compiles without complaint on the mac and runs in about 15 seconds but the check vector is always full of NaN's

  4. #4
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152

    Re: C --> Timing a Matrix Inversion

    [font=courier new][size=8pt]I don't know where the NaNs are coming from, jack. When I run it, the check vector is always filled with values, each of which is very very close to 1.

    From what you wrote, I am guessing that the check vector was correct when you ran it on your Mac using XP, is that true? I hope so.

    But, when you compiled it with a Mac C compiler, you got the bad check vectors? If so, I don't know why that would be.


    Dan

    (I see that wxDev-C++ is still being developed, Kent; --> that's good. But, it seems that nothing has been done with Dev-C++, since 2005; --> maybe, that's not so good. On the other hand, I guess C and C++ have not changed much lately.)
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  5. #5

    Re: C --> Timing a Matrix Inversion

    Dan, yes to both of your questions, I will have to take a look at code.

  6. #6

    Re: C --> Timing a Matrix Inversion

    Dan, the problem on the mac is that the random numbers were 32-bit instead of 16.

  7. #7
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152

    Re: C --> Timing a Matrix Inversion

    [font=courier new][size=8pt]I put in the code for the output file, that the random numbers vary from 0 to 32767. I should have added that this range was for my machine (I checked RAND_MAX for my machine.), and not necessarily for every machine. According to "C in a Nutshell", RAND_MAX is equal to at least 32767, but, it can be bigger. It varies according to the machine architecture, I guess. But, the matrix that the random numbers are assigned to, is of type double (64 bits). So, whether the random numbers are 16 bits or 32 bits, I think the compiler should automatically convert them to doubles. Unless, for your Mac C compiler, you have to provide an explicit cast in the code. In that case, I guess the line in the function initializeall(), would change from,

    m1{i}{j} = rand();

    to,

    m1{i}{j} = (double) rand();

    (I substituted brackets ({}) for braces (&#91;]), above. The output of the forum editor doesn't display braces properly, probably because they are used to format text.)

    Concerning the explicit cast, I could be wrong, I'm no expert.

    -------------------------------------------------

    Anyway, I'm glad that you at least got the program to work on your Mac, using the virtual XP environment.

    :P
    Dan
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  8. #8
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152

    Re: C --> Timing a Matrix Inversion

    [font=courier new][size=8pt]I was able to reduce the time it takes to fill a 1000 x 1000 matrix with random
    doubles, and to invert it. On my machine, the reduction is approximately 26%
    (previously, the operation averaged 35 seconds on my machine, and now, it
    averages 27 seconds.)

    :P
    Dan

    [code=c]//------------------------------------------------------------------------------------------------

    // A 'C' program to time matrix inversions.
    // Uses the Gauss-Jordan Method to invert an N x N matrix, filled with random doubles.
    // The results are written to a file; by default, "matinv.txt".

    //------------------------------------------------------------------------------------------------

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

    //------------------------------------------------------------------------------------------------

    // Set size of N x N matrix to be inverted.
    #define N 1000

    //------------------------------------------------------------------------------------------------

    const int first = 0;
    const int last = N - 1;

    // Name of outfile.
    const char outfilename[80] = "matinv.txt";

    //------------------------------------------------------------------------------------------------

    void invert(double original[N][N], double inversion[N][N]);
    void initializeall(double m1[N][N], double m2[N][N], double m3[N][N], double v1[N]);
    void matmatmult(double m1[N][N], double m2[N][N], double m3[N][N]);
    void matvecmult(double m1[N][N], double v1[N], double v2[N]);
    double getelapsedtime();
    void writeoutfile(double t, double v1[N]);

    //------------------------------------------------------------------------------------------------

    // Matrix to be inverted (filled with random doubles).
    double a0[N][N];

    // Copy of matrix to be inverted.
    double am[N][N];

    // Calculated inversion of a0.
    double ai[N][N];

    // Calculated identity matrix.
    double id[N][N];

    // va is multiplied by id, to give the check vector, vb.
    double va[N];

    // check vector.
    double vb[N];

    //------------------------------------------------------------------------------------------------

    int main(int argc, char *argv&#91;])
    {
    double t0, t1, tt;

    printf("Inverting matrix..\n");
    printf("\n");

    t0 = getelapsedtime();

    // a0 is filled with random doubles.
    // am is set equal to a0.
    // ai is set to the N x N identity matrix.
    // va is set to an N-dimensional vector, of all '1's.
    initializeall(a0, am, ai, va);

    // am becomes the N x N identity matrix.
    // ai becomes the (N x N) inversion of a0.
    invert(am, ai);

    t1 = getelapsedtime();
    tt = t1 - t0;

    printf("Finished inverting a %d by %d matrix.\n", N, N);
    printf("elapsed time = %6.3f seconds\n", tt);
    printf("\n");
    printf("Calculating check vector..\n");
    printf("\n");

    // Multiply the inverted matrix, ai, by the original matrix, a0, to give the calculated identity matrix, id.
    matmatmult(a0, ai, id);

    // Multiply the N-dimensional vector of all '1's, va, by the calculated N x N identity matrix, id, to give the check vector, vb.
    // If the inversion is correct, then, vb, should also be an N-dimensional vector, of all '1's.
    // vb, is written to the output file.
    matvecmult(id, va, vb);

    writeoutfile(tt, vb);

    printf("Finished calculating check vector.\n");
    printf("Output written to file, \"%s\".\n", outfilename);
    printf("\n");
    printf("Press 'Enter' to quit.\n");
    getchar();
    return 0;
    }

    //------------------------------------------------------------------------------------------------

    void invert(double original[N][N], double inversion[N][N])
    // Uses the Gauss-Jordan Method.
    {
    int pivot, maxrow, row, col;
    double max, test, temp, divisor, factor;
    //---------------------------------------------
    for(pivot = first; pivot <= last; pivot++)
    {
    //---------------------------------------------
    for(row = pivot; row <= last; row++)
    {
    max = 0;
    test = abs(original[row][pivot]);
    if( test > max)
    {
    max = test;
    maxrow = row;
    }
    }
    //---------------------------------------------
    if(maxrow != pivot)
    {
    for(col = first; col <= last; col++)
    {
    temp = original[pivot][col];
    original[pivot][col] = original[maxrow][col];
    original[maxrow][col] = temp;
    temp = inversion[pivot][col];
    inversion[pivot][col] = inversion[maxrow][col];
    inversion[maxrow][col] = temp;
    }
    }
    //---------------------------------------------
    divisor = original[pivot][pivot];
    for(col = first; col <= last; col++)
    {
    original[pivot][col] /= divisor;
    inversion[pivot][col] /= divisor;
    }
    //---------------------------------------------
    for(row = first; row <= last; row++)
    {
    if(row != pivot)
    {
    factor = -original[row][pivot];
    for(col = first; col <= last; col++)
    {
    original[row][col] += original[pivot][col] * factor;
    inversion[row][col] += inversion[pivot][col] * factor;
    }
    }
    }
    }
    }

    //------------------------------------------------------------------------------------------------

    void initializeall(double m1[N][N], double m2[N][N], double m3[N][N], double v1[N])
    {
    int i, j;
    time_t t;

    srand(time(&t));

    for(i = first; i <= last; i++)
    for(j = first; j <= last; j++)
    {
    m1[i][j] = rand();
    m2[i][j] = m1[i][j];
    }

    for(i = first; i <= last; i++)
    {
    v1[i] = 1;
    for(j = first; j <= last; j++)
    {
    m3[i][j] = 0;
    if(i == j) m3[i][j] = 1;
    }
    }
    }

    //------------------------------------------------------------------------------------------------

    void matmatmult(double m1[N][N], double m2[N][N], double m3[N][N])
    {
    int i, j, k;
    for(i = 0; i <= last; i++)
    {
    for(j = 0; j <= last; j++)
    {
    m3[i][j] = 0;
    for(k = 0; k <= last; k++) m3[i][j] += m1[i][k] * m2[k][j];
    }
    }
    }

    //------------------------------------------------------------------------------------------------

    void matvecmult(double m1[N][N], double v1[N], double v2[N])
    {
    int i, j;
    for(i = 0; i <= last; i++)
    {
    v2[i] = 0;
    for(j = 0; j <= last; j++) v2[i] += m1[i][j] * v1[j];
    }
    }

    //------------------------------------------------------------------------------------------------

    double getelapsedtime()
    {
    time_t ticks;
    double et;
    ticks = clock();
    et = (double)ticks / CLOCKS_PER_SEC;
    return et;
    }

    //------------------------------------------------------------------------------------------------

    void writeoutfile(double t, double v1[N])
    {
    int i, j;
    FILE *outfile;

    outfile = fopen(outfilename, "w");
    fprintf(outfile, "Inversion of a %d by %d matrix,\n", N, N);
    fprintf(outfile, "filled with uniformly distributed random doubles,\n");
    fprintf(outfile, "in the range of, [0, 32767] (on my machine).\n");
    fprintf(outfile, "\n");
    fprintf(outfile, "elapsed time = %6.3f seconds\n", t);
    fprintf(outfile, "\n");
    fprintf(outfile, "If the inversion was calculated correctly, then every\n");
    fprintf(outfile, "value shown below should be very close to 1.\n");
    fprintf(outfile, "\n");

    for(i = first; i <= last; i++) fprintf(outfile, "%5i %19.16f\n", i + 1, v1[i]);

    fclose(outfile);
    }

    //------------------------------------------------------------------------------------------------
    [/code]
    Attached Files Attached Files
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  9. #9
    thinBasic MVPs kryton9's Avatar
    Join Date
    Nov 2006
    Location
    Naples, Florida & Duluth, Georgia
    Age
    67
    Posts
    3,869
    Rep Power
    404

    Re: C --> Timing a Matrix Inversion

    You guys got some fast PC's. I get 30.14 secs for the first program and 29.56 for the second.
    Acer Notebook: Win 10 Home 64 Bit, Core i7-4702MQ @ 2.2Ghz, 12 GB RAM, nVidia GTX 760M and Intel HD 4600
    Raspberry Pi 3: Raspbian OS use for Home Samba Server and Test HTTP Server

  10. #10

    Re: C --> Timing a Matrix Inversion

    Dan, I compiled your first example with gcc with the command
    gcc main.c -o main
    execution time on the mac was about 15 seconds, with optimization it was about 7 seconds, under windows the speedup was only a few seconds.
    gcc main.c -O3 -o main

    under Windows XP(32) your second example runs in about 18 second unoptimized and about 13 seconds optimized
    on the Mac it runs in about 12 and 7 seconds repectivly.

Page 1 of 2 12 LastLast

Members who have read this thread: 0

There are no members to list at the moment.

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •