C++

Forum

 
 
How To Run SGESV Program?

Pages: 12
kindgnice (120)
Hello Profesionals,

Good day. I am trying to run this kind of program from Intel.

https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/sgesv_ex.c.htm

Sample Program:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
 
Intel® Math Kernel Library LAPACK Examples

/*******************************************************************************
*  Copyright (C) 2009-2015 Intel Corporation. All Rights Reserved.
*  The information and material ("Material") provided below is owned by Intel
*  Corporation or its suppliers or licensors, and title to such Material remains
*  with Intel Corporation or its suppliers or licensors. The Material contains
*  proprietary information of Intel or its suppliers and licensors. The Material
*  is protected by worldwide copyright laws and treaty provisions. No part of
*  the Material may be copied, reproduced, published, uploaded, posted,
*  transmitted, or distributed in any way without Intel's prior express written
*  permission. No license under any patent, copyright or other intellectual
*  property rights in the Material is granted to or conferred upon you, either
*  expressly, by implication, inducement, estoppel or otherwise. Any license
*  under such intellectual property rights must be express and approved by Intel
*  in writing.
*
********************************************************************************
*/
/*
   SGESV Example.
   ==============
 
   The program computes the solution to the system of linear
   equations with a square matrix A and multiple
   right-hand sides B, where A is the coefficient matrix:
 
     6.80  -6.05  -0.45   8.32  -9.67
    -2.11  -3.30   2.58   2.71  -5.14
     5.66   5.36  -2.70   4.35  -7.26
     5.97  -4.44   0.27  -7.17   6.08
     8.23   1.08   9.04   2.14  -6.87

   and B is the right-hand side matrix:
 
     4.02  -1.56   9.81
     6.19   4.00  -4.09
    -8.22  -8.67  -4.57
    -7.57   1.75  -8.61
    -3.03   2.86   8.99
 
   Description.
   ============
 
   The routine solves for X the system of linear equations A*X = B,
   where A is an n-by-n matrix, the columns of matrix B are individual
   right-hand sides, and the columns of X are the corresponding
   solutions.

   The LU decomposition with partial pivoting and row interchanges is
   used to factor A as A = P*L*U, where P is a permutation matrix, L
   is unit lower triangular, and U is upper triangular. The factored
   form of A is then used to solve the system of equations A*X = B.

   Example Program Results.
   ========================
 
 SGESV Example Program Results

 Solution
  -0.80  -0.39   0.96
  -0.70  -0.55   0.22
   0.59   0.84   1.90
   1.32  -0.10   5.36
   0.57   0.11   4.04

 Details of LU factorization
   8.23   1.08   9.04   2.14  -6.87
   0.83  -6.94  -7.92   6.55  -3.99
   0.69  -0.67 -14.18   7.24  -5.19
   0.73   0.75   0.02 -13.82  14.19
  -0.26   0.44  -0.59  -0.34  -3.43

 Pivot indices
      5      5      3      4      5
*/
#include <stdlib.h>
#include <stdio.h>

/* SGESV prototype */
extern void sgesv( int* n, int* nrhs, float* a, int* lda, int* ipiv,
                float* b, int* ldb, int* info );
/* Auxiliary routines prototypes */
extern void print_matrix( char* desc, int m, int n, float* a, int lda );
extern void print_int_vector( char* desc, int n, int* a );

/* Parameters */
#define N 5
#define NRHS 3
#define LDA N
#define LDB N

/* Main program */
int main() {
        /* Locals */
        int n = N, nrhs = NRHS, lda = LDA, ldb = LDB, info;
        /* Local arrays */
        int ipiv[N];
        float a[LDA*N] = {
            6.80f, -2.11f,  5.66f,  5.97f,  8.23f,
           -6.05f, -3.30f,  5.36f, -4.44f,  1.08f,
           -0.45f,  2.58f, -2.70f,  0.27f,  9.04f,
            8.32f,  2.71f,  4.35f, -7.17f,  2.14f,
           -9.67f, -5.14f, -7.26f,  6.08f, -6.87f
        };
        float b[LDB*NRHS] = {
            4.02f,  6.19f, -8.22f, -7.57f, -3.03f,
           -1.56f,  4.00f, -8.67f,  1.75f,  2.86f,
            9.81f, -4.09f, -4.57f, -8.61f,  8.99f
        };
        /* Executable statements */
        printf( " SGESV Example Program Results\n" );
        /* Solve the equations A*X = B */
        sgesv( &n, &nrhs, a, &lda, ipiv, b, &ldb, &info );
        /* Check for the exact singularity */
        if( info > 0 ) {
                printf( "The diagonal element of the triangular factor of A,\n" );
                printf( "U(%i,%i) is zero, so that A is singular;\n", info, info );
                printf( "the solution could not be computed.\n" );
                exit( 1 );
        }
        /* Print solution */
        print_matrix( "Solution", n, nrhs, b, ldb );
        /* Print details of LU factorization */
        print_matrix( "Details of LU factorization", n, n, a, lda );
        /* Print pivot indices */
        print_int_vector( "Pivot indices", n, ipiv );
        exit( 0 );
} /* End of SGESV Example */

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, float* a, int lda ) {
        int i, j;
        printf( "\n %s\n", desc );
        for( i = 0; i < m; i++ ) {
                for( j = 0; j < n; j++ ) printf( " %6.2f", a[i+j*lda] );
                printf( "\n" );
        }
}

/* Auxiliary routine: printing a vector of integers */
void print_int_vector( char* desc, int n, int* a ) {
        int j;
        printf( "\n %s\n", desc );
        for( j = 0; j < n; j++ ) printf( " %6i", a[j] );
        printf( "\n" );
}


However, when I tried to run it using Microsoft Visual Studio, it says:

Link: https://i.imgur.com/oMKvmVQ.png

What can I do in order to make it work?

Thank you for your inputs in advance. God bless.
Last edited on
lastchance (6980)
Looks like you haven't linked the library that has routine sgesv() in.

Might be quicker to write your own LU factoriser.
Last edited on
kindgnice (120)
Thanks for your response @lastchance. Oh, do you mean it's "more faster" to use the typical linear solver rather than using this library?
lastchance (6980)
@kindgnice,
To use this routine as is, you will need to specify whic library to use (consult the documentation for this routine) and let the Visual Studio IDE know that this has to be linked (for which you will need Visual Studio documentation). I don't use MS Visual Studio (or LaPack), so I can't advise you on either.

I have to ask, in the light of your other current threads, WHAT are you actually intending to do with this routine? Maybe you don't need it. It is relatively straightforward to do an LU factorisation for small, amenable matrices, but less so if it needs partial pivoting or for large matrices.
kindgnice (120)
Thank you for your response @lastchance. You are correct @lastchance. After I downloaded the "clapack.h", and tried running again the program, it worked. :)

I am still new with this SGESV function, but what I have heard from my professor is that, this SGESV function can be used in order to solve linear equations such as Ax=B, wherein x is unknown. But still, I am still new to this and trying to digest it up to now... :(

Some of my classmates also told me that this SGESV function has some untypical behavior due to the fact that the display of matrix from clapack is kinda different because it was built in fortran, but this information is something that I have to figure out as well... :(
Last edited on
lastchance (6980)
Yes, the conversion from 1d array to i,j positions implied by line 136 is the Fortran ordering convention: go down first column, then second column, etc. C and related languages go across first row, then second row etc.
kindgnice (120)
Thank you for expounding the details @lastchance.

Unfortunately, I thought I already successfully run the program by doing this:


1
2
3
4
5
6
7
8
9
10
11
12
13
 
#include <iostream>
#include <fstream>
#include <cstdlib> // exit()
#include <cmath> // floorf()

extern "C" {
#include "f2c.h"
#include "clapack.h"
}

int sgesv_(integer * n, integer * nrhs, real * a, integer * lda, integer * ipiv, real * b, integer * ldb, integer * info) {
	return 0;
}


I thought initializing sgesv_ by creating a function will do the trick. It's running, yes, until I finally realized that sgesv_ is not totally working for me because I am just returning 0 values?

:(

P.S. I guess I still have no access to sgesv_ function yet even though I already included the "clapack" library header.
Last edited on
lastchance (6980)
You don't need to "initialise" sgesv. It's just a function in an external library. Just call it when needed.

As far as I can see from the example given in the initial post, sgesv (WITHOUT THE UNDERSCORE _) is the routine in the library. Haven't you run that example?

And we still come back to the question: what are you intending to do with your matrix solver at the end of the day? It doesn't seem to have much application in your other threads.


If you want a cheap LU factoriser here's one. NOTE: it WON'T WORK FOR ALL MATRICES - it doesn't have the partial pivoting of the more advanced LAPACK solver, and it probably won't be stable for large matrices.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
 
//======================================================================
// Here, [m x n] matrices are represented as 1-d arrays
// Row i goes from 0 to m-1.
// Col j goes from 0 to n-1.
// The (i,j) component is stored at 1-d index
//                 n * i + j
//
// The storage order is the normal one for C++ arrays (column index varies fastest);
//   ie across top of the matrix first
// NOTE: opposite storage order to Fortran
//======================================================================

#include <iostream>
#include <sstream>
#include <iomanip>
#include <cmath>
using namespace std;

//======================================================================

template <class T> void LUFactorise( int n, T *A, T *L, T *U )       // Factorise A = LU, where L and U are
{                                                                    //    lower- and upper-triangular respectively
   // Initialise
   for ( int ij = 0; ij < n * n; ij++ ) L[ij] = U[ij] = 0.0;

   for ( int i = 0; i < n; i++ )
   {
      // ith row -> U
      for ( int j = i; j < n; j++ )
      {
         int ij = n * i + j;
         U[ij] = A[ij];
         for ( int k = 0; k < i; k++ ) U[ij] -= L[n*i+k] * U[n*k+j];
      }

      // ith column -> L
      for ( int j = i; j < n; j++ )
      {
         int ji = n * j + i;
         L[ji] = A[ji];
         for ( int k = 0; k < i; k++ ) L[ji] -= L[n*j+k] * U[n*k+i];
         L[ji] /= U[n*i+i];
      }
   }
}

//======================================================================

template <class T> void solveUpper( int n, T *U, T *B, T *X )        // Solve UX = B, where U is upper-triangular
{
   for ( int i = n - 1; i >= 0; i-- )
   {
      X[i] = B[i];
      for ( int j = i + 1; j < n; j++ )  X[i] -= U[n*i+j] * X[j];
      X[i] /= U[n*i+i];
   }
}

//======================================================================

template <class T> void solveLower( int n, T *L, T *B, T *X )        // Solve LX = B, where L is lower-triangular
{
   for ( int i = 0; i < n; i++ )
   {
      X[i] = B[i];
      for ( int j = 0; j < i; j++ ) X[i] -= L[n*i+j] * X[j];
      X[i] /= L[n*i+i];
   }
}

//======================================================================

template <class T> void matMul( int mA, int nA, int nB, T *A, T *B, T *C )     // Matrix multiply: A * B = C
{
   for ( int i = 0; i < mA; i++ )
   {
      for ( int j = 0; j < nB; j++ )
      {
         int ij = nB * i + j;
         C[ij] = 0.0;
         for ( int k = 0; k < nA; k++ ) C[ij] += A[nA*i+k] * B[nB*k+j];
      }
   }
}

//======================================================================

template <class T> void showMatrix( int m, int n, T *A )
{
   const double NEARZERO = 1.0e-10;              // You will have to tailor
   const string SPACE = "  ";                    // all these formatting parameters
   const int W = 12;                             //
   const int P = 4;                              //
   cout << fixed << setprecision( P );
// cout << scientific << setprecision( P );

   for ( int i = 0; i < m; i++ )
   {
      for ( int j = 0; j < n; j++ )
      {
         T val = A[n*i+j];   if ( abs( val ) < NEARZERO ) val = 0.0;
         cout << setw( W ) << val << SPACE;
      }
      cout << '\n';
   }
}

//======================================================================

int main()
{
   int N = 5;
   double *A = new double[N*N];        // Matrix
   double *B = new double[N*1];        // RHS vector
   double *X = new double[N*1];        // Solution

   double *L = new double[N*N];        // Lower-triangular matrix
   double *U = new double[N*N];        // Upper-triangular matrix
   double *Y = new double[N*1];        // Temporary solution of LY = B  (whence UX = Y gives full solution)

   double *AX= new double[N*1];        // Check matrix AX (should equal B)


   stringstream inA( "  6.80 -2.11  5.66  5.97  8.23 "
                     " -6.05 -3.30  5.36 -4.44  1.08 "
                     " -0.45  2.58 -2.70  0.27  9.04 "
                     "  8.32  2.71  4.35 -7.17  2.14 "
                     " -9.67 -5.14 -7.26  6.08 -6.87 " );
   stringstream inB( "  4.02  6.19 -8.22 -7.57 -3.03 " );


   // Input data (here, just from stringstreams)
   for ( int ij = 0; ij < N * N; ij++ ) inA >> A[ij];
   for ( int i  = 0; i  < N    ; i++  ) inB >> B[i ];


   // LU factorise; then solve
   LUFactorise( N, A, L, U );
   solveLower( N, L, B, Y );
   solveUpper( N, U, Y, X );           // Final solution is in X 


   // Write summary
   cout << "\n*** Solve AX = B by LU decomposition ***\n";
   cout << "\nOriginal matrix A:        \n";    showMatrix( N, N, A );
   cout << "\nLower-triangular matrix L:\n";    showMatrix( N, N, L );
   cout << "\nUpper-triangular matrix U:\n";    showMatrix( N, N, U );
   cout << "\nSolution vector X:        \n";    showMatrix( N, 1, X );


   // Check solution
   matMul( N, N, 1, A, X, AX );
   cout << "\n\n*** CHECK ***\n";
   cout << "\nAX:\n";    showMatrix( N, 1, AX );
   cout << "\nB :\n";    showMatrix( N, 1, B  );


   delete [] A;
   delete [] B;
   delete [] X;
   delete [] L;
   delete [] U;
   delete [] Y;
   delete [] AX;
}



*** Solve AX = B by LU decomposition ***

Original matrix A:        
      6.8000       -2.1100        5.6600        5.9700        8.2300  
     -6.0500       -3.3000        5.3600       -4.4400        1.0800  
     -0.4500        2.5800       -2.7000        0.2700        9.0400  
      8.3200        2.7100        4.3500       -7.1700        2.1400  
     -9.6700       -5.1400       -7.2600        6.0800       -6.8700  

Lower-triangular matrix L:
      1.0000        0.0000        0.0000        0.0000        0.0000  
     -0.8897        1.0000        0.0000        0.0000        0.0000  
     -0.0662       -0.4714        1.0000        0.0000        0.0000  
      1.2235       -1.0221        3.1267        1.0000        0.0000  
     -1.4221        1.5724       -6.0422       -1.1624        1.0000  

Upper-triangular matrix U:
      6.8000       -2.1100        5.6600        5.9700        8.2300  
      0.0000       -5.1773       10.3957        0.8715        8.4023  
      0.0000        0.0000        2.5747        1.0759       13.5451  
      0.0000        0.0000        0.0000      -16.9476      -41.6927  
      0.0000        0.0000        0.0000        0.0000       25.0011  

Solution vector X:        
     -1.4658  
      2.7127  
      2.8994  
      1.8569  
     -0.9460  


*** CHECK ***

AX:
      4.0200  
      6.1900  
     -8.2200  
     -7.5700  
     -3.0300  

B :
      4.0200  
      6.1900  
     -8.2200  
     -7.5700  
     -3.0300

Last edited on
kindgnice (120)
Thank you very much for your excellent support @lastchance, as always. Apologies if it took me so long to respond. I was a bit pre-occupied of other subjects due to presentations and stuffs. I found some articles about Crout, Doolittle and Cholesky factorization. I wonder if those are efficient or better than SGESV?
Last edited on
kindgnice (120)
Apparently, when I tried to use the code above, I got these results:

Given:



stringstream inA(" 0 0.3019 0.1562 0.027 1 0.2 0.5 "
" 0.3019 0 0.1562 0.7637 1 0.8 0.2 "
" 0.1562 0.1326 0 0.02263 1 0.7 0.7 "
" 0.027 0.7637 0.02263 0 1 0.5 0.5 "
" 1 1 1 1 0 0 0 "
" 0.2 0.8 0.7 0.5 0 0 0 "
" 0.5 0.2 0.7 0.5 0 0 0 "
);
stringstream inB(" 0.0 0.0 0.0 0.2 0.0 0.0 0.0 ");



My output:


Original matrix A:
      0.0000        0.3019        0.1562        0.0270        1.0000        0.2000        0.5000
      0.3019        0.0000        0.1562        0.7637        1.0000        0.8000        0.2000
      0.1562        0.1326        0.0000        0.0226        1.0000        0.7000        0.7000
      0.0270        0.7637        0.0226        0.0000        1.0000        0.5000        0.5000
      1.0000        1.0000        1.0000        1.0000        0.0000        0.0000        0.0000
      0.2000        0.8000        0.7000        0.5000        0.0000        0.0000        0.0000
      0.5000        0.2000        0.7000        0.5000        0.0000        0.0000        0.0000

Lower-triangular matrix L:
   -nan(ind)        0.0000        0.0000        0.0000        0.0000        0.0000        0.0000
         inf     -nan(ind)        0.0000        0.0000        0.0000        0.0000        0.0000
         inf     -nan(ind)     -nan(ind)        0.0000        0.0000        0.0000        0.0000
         inf     -nan(ind)     -nan(ind)     -nan(ind)        0.0000        0.0000        0.0000
         inf     -nan(ind)     -nan(ind)     -nan(ind)     -nan(ind)        0.0000        0.0000
         inf     -nan(ind)     -nan(ind)     -nan(ind)     -nan(ind)     -nan(ind)        0.0000
         inf     -nan(ind)     -nan(ind)     -nan(ind)     -nan(ind)     -nan(ind)     -nan(ind)

Upper-triangular matrix U:
      0.0000        0.3019        0.1562        0.0270        1.0000        0.2000        0.5000
      0.0000          -inf          -inf          -inf          -inf          -inf          -inf
      0.0000        0.0000     -nan(ind)     -nan(ind)     -nan(ind)     -nan(ind)     -nan(ind)
      0.0000        0.0000        0.0000     -nan(ind)     -nan(ind)     -nan(ind)     -nan(ind)
      0.0000        0.0000        0.0000        0.0000     -nan(ind)     -nan(ind)     -nan(ind)
      0.0000        0.0000        0.0000        0.0000        0.0000     -nan(ind)     -nan(ind)
      0.0000        0.0000        0.0000        0.0000        0.0000        0.0000     -nan(ind)

Solution vector X:
   -nan(ind)
   -nan(ind)
   -nan(ind)
   -nan(ind)
   -nan(ind)
   -nan(ind)
   -nan(ind)


*** CHECK ***

AX:
   -nan(ind)
   -nan(ind)
   -nan(ind)
   -nan(ind)
   -nan(ind)
   -nan(ind)
   -nan(ind)

B :
      0.0000
      0.0000
      0.0000
      0.2000
      0.0000
      0.0000
      0.0000
kindgnice (120)
I don't know why the answers are not solved...

However, if I use this linear solver that was shared to me by my professor, at some point, it works for any values though...

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
 
	// ==============================================
	// == Start linear solver (LU decomposition)
	// ==============================================
	int p[NMAX], r = 0;
	float amax, w, akk1, ptemp, atemp;

	for (int i = 0; i < n; i++) /* Index for Pivoting */
		p[i] = i;
	/* Forward Substitution with Pivoting */
	for (int k = 0; k < n - 1; k++) {
		amax = fabs(aLU[k][k]);
		r = k;
		for (int i = k + 1; i < n; i++) { /* search r such that a_{r,k}=max a_{i,k} i=k,..n */
			if (fabs(aLU[i][k]) > amax) {
				amax = fabs(aLU[i][k]);
				r = i;
			}
		}
		if (r != k) { /* exchange column r and column k */
			ptemp = p[r];
			p[r] = p[k];
			p[k] = ptemp;
			for (int j = 0; j < n + 1; j++) {
				atemp = aLU[r][j];
				aLU[r][j] = aLU[k][j];
				aLU[k][j] = atemp;
			}
		}
		for (int i = k + 1; i < n; i++)
			aLU[i][k] = aLU[i][k] / aLU[k][k];
		for (int i = k + 1; i < n; i++) {
			for (int j = k + 1; j < n + 1; j++)
				aLU[i][j] -= aLU[i][k] * aLU[k][j];
		}
	}
	/* Backward Substitution */
	for (int k = n - 1; k >= 0; --k) {
		akk1 = 1 / aLU[k][k];
		for (int i = k + 1; i < n; i++) {
			aLU[k][n] -= aLU[k][i] * aLU[i][n]; //before: aLU[k][n] -= aLU[k][i] * aLU[i][n];
		}
		aLU[k][n] *= akk1;
	}
	// ==============================================
	// == End linear solver (LU decomposition)
	// ==============================================
kindgnice (120)
For your question why I need LU decomposition...

I need this because there are lots of times I encounter matrixes in the form of Ax=B, and this LU decomposition codes will help me find the x values.
lastchance (6980)
Yes, well that's why SGESV includes partial pivoting as well as LU factorisation. You have supplied a matrix with 0 in the [0][0] position, so it won't be able to use the first pivot without swapping rows.

Swap the first two rows of A (and B, although that won't matter in your example). The equation set will still be the same, just in a different order. Then it will work. This swapping of rows is what partial pivoting does (automatically); I haven't implemented it in my code yet, as I have no particular need for such a solver at present.

EDIT: our answers crossed; your professor's code includes partial pivoting (searching out the largest element on or below the pivot and then swapping rows(?)). You can use that.



In reply to one of your previous posts you aren't comparing like with like. LU, Crout's, Doolittle, Cholesky's method are ALGORITHMS for solving problems; SGESV is a LIBRARY routine that happens to implement one of them.

Which method or library routine depends on your problem. Do you want a single matrix equation solved (in which case better to use a direct method) or with several right-hand sides (in which case use a factorisation method or even invert the matrix). Is the matrix large or small? Does the matrix have a special form? (e.g. tridiagonal; symmetric). Is the matrix sparse? Without saying what you intend to do it would be impossible to answer your question.

Note that Cholesky factorisation only works for symmetric and positive-definite matrices (or Hermitian matrices if they are complex rather than real). If your matrix is symmetric and positive definite then it is better than LU factorisation; if not, then it won't work. So you need to consider your actual application.

Unless you want to write your own solvers I think you should spend time making sure that you can run whatever linear-algebra library routines you have available and note the form in which you supply matrices and vectors to them.



For the record, here's what you get if you first swap the first two rows of A and the first two rows of B before applying factorisation. Obviously, I have done this swap manually: partial pivoting would do it automatically. It would also be more stable for large matrices. I haven't time to implement the necessary row-swapping at present. If your professor's code works then just use that.


EDITED - forgot to update N
*** Solve AX = B by LU decomposition ***

Original matrix A:        
      0.3019        0.0000        0.1562        0.7637        1.0000        0.8000        0.2000  
      0.0000        0.3019        0.1562        0.0270        1.0000        0.2000        0.5000  
      0.1562        0.1326        0.0000        0.0226        1.0000        0.7000        0.7000  
      0.0270        0.7637        0.0226        0.0000        1.0000        0.5000        0.5000  
      1.0000        1.0000        1.0000        1.0000        0.0000        0.0000        0.0000  
      0.2000        0.8000        0.7000        0.5000        0.0000        0.0000        0.0000  
      0.5000        0.2000        0.7000        0.5000        0.0000        0.0000        0.0000  

Lower-triangular matrix L:
      1.0000        0.0000        0.0000        0.0000        0.0000        0.0000        0.0000  
      0.0000        1.0000        0.0000        0.0000        0.0000        0.0000        0.0000  
      0.5174        0.4392        1.0000        0.0000        0.0000        0.0000        0.0000  
      0.0894        2.5296        2.5864        1.0000        0.0000        0.0000        0.0000  
      3.3124        3.3124        0.2328       -1.7838        1.0000        0.0000        0.0000  
      0.6625        2.6499       -1.2221       -0.6381        0.4488        1.0000        0.0000  
      1.6562        0.6625       -2.2609       -1.9261        0.5714        0.4765        1.0000  

Upper-triangular matrix U:
      0.3019        0.0000        0.1562        0.7637        1.0000        0.8000        0.2000  
      0.0000        0.3019        0.1562        0.0270        1.0000        0.2000        0.5000  
      0.0000        0.0000       -0.1494       -0.3844        0.0434        0.1982        0.3769  
      0.0000        0.0000        0.0000        0.8575       -1.7313       -0.5902       -1.7576  
      0.0000        0.0000        0.0000        0.0000       -9.7231       -4.4113       -5.5415  
      0.0000        0.0000        0.0000        0.0000        0.0000        0.7857        0.3689  
      0.0000        0.0000        0.0000        0.0000        0.0000        0.0000       -0.2053  

Solution vector X:        
      0.3539  
      0.1769  
      0.2654  
     -0.7962  
      0.3244  
      0.4090  
     -0.9592  


*** CHECK ***

AX:
      0.0000  
      0.0000  
      0.0000  
      0.2000  
      0.0000  
      0.0000  
      0.0000  

B :
      0.0000  
      0.0000  
      0.0000  
      0.2000  
      0.0000  
      0.0000  
      0.0000
Last edited on
kindgnice (120)
Thank you for the swift response @lastchance. I think I could consider my matrices as sparse because for every run of my program, there are lots of zeros in it. As I will be working with more complex data points (such as bunny, armadillo, in the future for example), I think I need a more flexible LU decomposition solver for bigger matrices that can handle the coordinates of bunny, armadillo, etc.

Taking into consideration the generated matrices I had so far, especially for my Matrix A (of Ax=B equations), it looks symmetric. The matrix of A has always nxn rows and columns.
kindgnice (120)
In addition to the program I shared above, if I will try to run it, the output that I will get will be:


Matrix A =
0       0.3019  0.1562  0.027   1       0.2     0.5
0.3019  0       0.1326  0.07637 1       0.8     0.2
0.1562  0.1326  0       0.02263 1       0.7     0.7
0.027   0.07637 0.02263 0       1       0.5     0.5
1       1       1       1       0       0       0
0.2     0.8     0.7     0.5     0       0       0
0.5     0.2     0.7     0.5     0       0       0

Vector b =
0
0
0
0.2
0
0
0

Solution x =
-1.693
-0.8463
-1.269
3.808
0.4608
-0.03948
-0.204
lastchance (6980)
kindgnice wrote:
it looks symmetric. The matrix of A has always nxn rows and columns.

That's not what makes a matrix symmetric, @kindgnice. "Symmetric" means that the matrix is the same if reflected in the main diagonal; i.e. it is equal to its transpose.

If your professor has shared working code with you then I should use that: it actually looks more usable than Armadillo. Your professor would also more easily be able to give you advice about his own code which will make your task much easier.
kindgnice (120)
I forgot to mention, this code that he shared to me was just taken from the internet too. The reason why I am trying to look for other LU method because I am trying to make a program in linear arrays way. The code above that the professor got from internet uses two dimensional arrays. And it has some problems especially for larger matrices. In this reason, I revamped my code (entirely) and transformed it into linear arrays (i.e. array[i*n+j]).

I think using linear arrays has less cost compared to using two dimensional arrays. Also, two dimensional LU decomposition versions have difficulty in large matrices.

And, I think you are right, the LU decomposition should also have a pivoting feature... Unfortunately, the codes from the link I saw using Cholesky, Crout and Doolittle has no pivoting I think...

http://www.sci.utah.edu/~wallstedt/LU.htm
Last edited on
lastchance (6980)
I'll have a look at the partial pivoting and see if I can replicate what SGESV would do. However, this is unlikely to be in the next day or two.

The partial pivoting is supposed to alleviate the difficulties with large matrices.

If your matrices are tridiagonal then there are definitely faster algorithms ...
kindgnice (120)
Thank you very much Sir in advance.
I am so sorry that I am so much incompetent to resolve this with myself. :(
Thank you very much and God bless....
lastchance (6980)
Try this please, @kindgnice.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
 
//======================================
// [m x n] matrices are represented as 1-d arrays
// Row i goes from 0 to m-1.
// Col j goes from 0 to n-1.
// The (i,j) component is stored at 1-d index
//    n * i + j
//
// Storage order as for C++ arrays (column index varies fastest);
// NOTE: opposite storage order to Fortran
//======================================
// Factorise A = PLU, where:
//    L and U are lower and upper triangular matrices
//    P is a permutation matrix
// P is returned here as an n-element array,
//    such that P[i] is the row that must be swapped into the ith row
// When solving Ax=b, do this swapping to b first, then solve LUx=b'
//======================================

#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <cmath>
#include <algorithm>
using namespace std;

const bool doPivot = true;
const double EPS = 1.0e-30;

//======================================


template <class T> bool LUFactorise( int n, T *A, T *L, T *U, int *P )
{                                                                   
   T *AA = new T[n*n];

   // Initialise
   fill( L, L + n * n, 0  );
   fill( U, U + n * n, 0  );
   copy( A, A + n * n, AA );
   for ( int i = 0; i < n; i++ ) P[i] = i;


   for ( int i = 0; i < n; i++ )
   {
      int r = i;                              // Find row r below with largest element in column i
      T maxA = abs( AA[n*i+i] );
      for ( int k = i + 1; k < n; k++ )       
      {
         int ki = n * k + i;
         if ( abs( AA[ki] ) > maxA )
         {
            r = k;
            maxA = abs( AA[ki] );
         }
      }

      if ( doPivot && r != i )
      {
         swap( P[i], P[r] );
         for ( int j = 0; j < n; j++ ) 
         {
            int rj = n * r + j;
            int ij = n * i + j;
            swap( AA[rj], AA[ij] );
            if ( j < i ) swap( L[rj], L[ij] );
         }
      }

      // ith row -> U
      for ( int j = i; j < n; j++ )
      {
         int ij = n * i + j;
         U[ij] = AA[ij];
         for ( int k = 0; k < i; k++ ) U[ij] -= L[n*i+k] * U[n*k+j];
      }

      if ( abs( U[n*i+i] ) < EPS ) return false;

      // ith column -> L
      for ( int j = i; j < n; j++ )
      {
         int ji = n * j + i;
         L[ji] = AA[ji];
         for ( int k = 0; k < i; k++ ) L[ji] -= L[n*j+k] * U[n*k+i];
         L[ji] /= U[n*i+i];
      }
   }

   delete [] AA;
   return true;
}

//======================================

template <class T> void solveUpper( int n, T *U, T *B, T *X )        // Solve UX = B, where U is upper-triangular
{
   for ( int i = n - 1; i >= 0; i-- )
   {
      X[i] = B[i];
      for ( int j = i + 1; j < n; j++ )  X[i] -= U[n*i+j] * X[j];
      X[i] /= U[n*i+i];
   }
}

//======================================

template <class T> void solveLower( int n, T *L, T *B, T *X )        // Solve LX = B, where L is lower-triangular
{
   for ( int i = 0; i < n; i++ )
   {
      X[i] = B[i];
      for ( int j = 0; j < i; j++ ) X[i] -= L[n*i+j] * X[j];
      X[i] /= L[n*i+i];
   }
}

//======================================

template <class T> void permute( int m, int n, int *P, T *A )        // Permute rows of an mxn matrix A
{
   T *AA = new T[m*n];
   copy( A, A + m * n, AA );

   for ( int i = 0; i < m; i++ )
   {
      int r = P[i];
      for ( int j = 0; j < n; j++ ) A[n*i+j] = AA[n*r+j];
   }

   delete [] AA;
}

//======================================

template <class T> void solvePLU( int n, int *P, T *L, T *U, T *B, T *X )      // Solve PLUx = b
{
   double *PB = new double[n*1];
   double *Y  = new double[n*1];       // Temporary solution of LY = B  (whence UX = Y gives full solution)
   copy( B, B + n, PB );
   permute( n, 1, P, PB );             // Permuted rows of B
   solveLower( n, L, PB, Y );
   solveUpper( n, U, Y , X );          // Final solution is in X
   delete [] PB;
   delete [] Y;
}

//======================================

template <class T> void matMul( int mA, int nA, int nB, T *A, T *B, T *C )     // Matrix multiply: A * B = C
{
   for ( int i = 0; i < mA; i++ )
   {
      for ( int j = 0; j < nB; j++ )
      {
         int ij = nB * i + j;
         C[ij] = 0.0;
         for ( int k = 0; k < nA; k++ ) C[ij] += A[nA*i+k] * B[nB*k+j];
      }
   }
}

//======================================

template <class T> void pMatrix( int n, int *P, T *A )               // Create a permutation matrix A from P
{
   fill( A, A + n * n, 0 );
   for ( int i = 0; i < n; i++ ) A[n*P[i]+i] = 1;
}

//======================================

template <class T> void showMatrix( int m, int n, T *A )
{
   const double ZERO = 1.0e-10;
   const string SPACE = "  ";
   const int w = 12;
   const int p = 4;
   cout << fixed << setprecision( p );

   for ( int i = 0; i < m; i++ )
   {
      for ( int j = 0; j < n; j++ )
      {
         T val = A[n*i+j];   if ( abs( val ) < ZERO ) val = 0.0;
         cout << setw( w ) << val << SPACE;
      }
      cout << '\n';
   }
}

//======================================

template <class T> void getData( int &n, T *&A, T *&B )
{
// ifstream in( "matrix.txt" );
   stringstream in(
      " 7 "                                         // N
      " 0      0.3019 0.1562  0.027   1 0.2 0.5 "   // <---A
      " 0.3019 0      0.1562  0.7637  1 0.8 0.2 "   //   .
      " 0.1562 0.1326 0       0.02263 1 0.7 0.7 "   //   .
      " 0.027  0.7637 0.02263 0       1 0.5 0.5 "   //   .
      " 1      1      1       1       0 0   0   "   //   .
      " 0.2    0.8    0.7     0.5     0 0   0   "   //   .
      " 0.5    0.2    0.7     0.5     0 0   0   "   //   .
      " 0.0    0.0    0.0     0.2     0 0   0   "   // <---B
   );
   in >> n;
   A = new double[n*n];
   B = new double[n*1];
   for ( int ij = 0; ij < n * n; ij++ ) in >> A[ij];
   for ( int i  = 0; i  < n    ; i++  ) in >> B[i ];
}

//======================================

int main()
{
   int N;
   double *A;                     // Matrix
   double *B;                     // RHS vector


   // INPUT
   getData( N, A, B );
   int    *P = new int   [N*1];   // Permutation of rows
   double *L = new double[N*N];   // Lower-triangular matrix
   double *U = new double[N*N];   // Upper-triangular matrix
   double *X = new double[N*1];   // Solution of Ax = b


   // FACTORISE
   bool ok = LUFactorise( N, A, L, U, P );
   if ( !ok ) { cout << "Can't factorise\n";   return 1; }


   // SOLVE
   solvePLU( N, P, L, U, B, X );


   // SUMMARY
   double *PM = new double[N*N];        // Permutation matrix
   pMatrix( N, P, PM );
   cout << "\n*** Solve AX = B by LU decomposition ***\n";
   cout << "\nOriginal A:        \n";    showMatrix( N, N, A  );
   cout << "\nPermutation PM:    \n";    showMatrix( N, N, PM );
   cout << "\nLower-triangular L:\n";    showMatrix( N, N, L  );
   cout << "\nUpper-triangular U:\n";    showMatrix( N, N, U  );
   cout << "\nSolution X:        \n";    showMatrix( N, 1, X  );


   // CHECKS
   double *PL = new double[N*N];
   double *PLU= new double[N*N];
   matMul( N, N, N, PM, L, PL  );
   matMul( N, N, N, PL, U, PLU );
   cout << "\n\n*** CHECK FACTORISATION, A = PLU ***\n";
   cout << "\nA:\n";    showMatrix( N, N, A   );
   cout << "\nPLU:\n";  showMatrix( N, N, PLU );
   delete [] PL;
   delete [] PLU;
   delete [] PM;

   double *AX = new double[N*1];
   matMul( N, N, 1, A, X, AX );
   cout << "\n\n*** CHECK SOLUTION, Ax = b ***\n";
   cout << "\nAX:\n";   showMatrix( N, 1, AX );
   cout << "\nB :\n";   showMatrix( N, 1, B  );
   delete [] AX;


   delete [] A;
   delete [] B;
   delete [] X;
   delete [] P;
   delete [] L;
   delete [] U;
}

Last edited on
Pages: 12