# include # include # include # include void dgeev_ ( char *jobvl, char *jobvr, int *n, double *A, int *lda, double *wr, double *wi, double *vl, int *ldvl, double *vr, int *ldvr, double *work, int *lwork, int *info ); void dgesv_ ( int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *lbd, int *info ); int main ( ); void dgeev_test ( ); void dgesv_test ( ); double *clement1_matrix ( int n ); void timestamp ( ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: lapack_test() tests lapack(). Discussion: This example shows how a C++ calling program can interact with the lapack() library, which is written in Fortran. Licensing: This code is distributed under the MIT license. Modified: 06 June 2024 Author: John Burkardt */ { timestamp ( ); printf ( "\n" ); printf ( "lapack_test():\n" ); printf ( " C version\n" ); printf ( " lapack() is a Fortran library of linear algebra utilities.\n" ); printf ( " A C++ program can call it, if careful.\n" ); dgeev_test ( ); dgesv_test ( ); /* Terminate. */ printf ( "\n" ); printf ( "lapack_test():\n" ); printf ( " Normal end of execution,\n" ); timestamp ( ); return 0; } /******************************************************************************/ void dgeev_test ( ) /******************************************************************************/ /* Purpose: dgeev_test() tests dgeev(). Discussion: dgeev() computes the eigenvalues and eigenvectors of a matrix A, where the matrix A is stored in general format. Licensing: This code is distributed under the MIT license. Modified: 06 June 2024 Author: John Burkardt */ { int n = 5; double *A; int i; int info; char jobvl; char jobvr; int lda = n; int ldvl = 1; int ldvr = 1; int lwork; double *vl; double *vr; double *wi; double *work; double *wr; printf ( "\n" ); printf ( "dgeev_test():\n" ); printf ( " dgeev() computes eigenvalues and eigenvectors of a matrix A,\n" ); printf ( " with matrix A in general storage.\n" ); A = clement1_matrix ( n ); jobvl = 'N'; jobvr = 'N'; lda = n; wr = ( double * ) malloc ( n * sizeof ( double ) ); wi = ( double * ) malloc ( n * sizeof ( double ) ); vl = NULL; ldvl = 1; vr = NULL; ldvr = 1; work = ( double * ) malloc ( 3 * n * sizeof ( double ) ); lwork = 3 * n; dgeev_ ( &jobvl, &jobvr, &n, A, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, &info ); printf ( "\n" ); printf ( " Computed eigenvalues lambda:\n" ); for ( i = 0; i < n; i++ ) { printf ( " %g + %g i\n", wr[i], wi[i] ); } free ( A ); free ( wi ); free ( work ); free ( wr ); return; } /******************************************************************************/ void dgesv_test ( ) /******************************************************************************/ /* Purpose: dgesv_test() tests dgesv(). Discussion: dgesv() computes the solution x of a linear system A*x=b, where the matrix A is stored in general format. Licensing: This code is distributed under the MIT license. Modified: 06 June 2024 Author: John Burkardt */ { int n = 3; double *A; double *b; int i; int info; int *ipvt; int lda = n; int ldb = n; int nrhs = 1; /* Matrix is given by columns. */ A = ( double * ) malloc ( n * n * sizeof ( double ) ); A[0+0*3] = 1.0; A[1+0*3] = 4.0; A[2+0*3] = 7.0; A[0+1*3] = 2.0; A[1+1*3] = 5.0; A[2+1*3] = 8.0; A[0+2*3] = 3.0; A[1+2*3] = 6.0; A[2+2*3] = 0.0; b = ( double * ) malloc ( n * sizeof ( double ) ); b[0] = 14.0; b[1] = 32.0; b[2] = 23.0; ipvt = ( int * ) malloc ( n * sizeof ( int ) ); printf ( "\n" ); printf ( "dgesv_test():\n" ); printf ( " dgesv() solves a linear system A*x=b,\n" ); printf ( " with matrix A in general storage.\n" ); dgesv_ ( &n, &nrhs, A, &lda, ipvt, b, &ldb, &info ); printf ( "\n" ); printf ( " Computed solution x:\n" ); for ( i = 0; i < n; i++ ) { printf ( " %g", b[i] ); } printf ( "\n" ); free ( A ); free ( b ); free ( ipvt ); return; } /******************************************************************************/ double *clement1_matrix ( int n ) /******************************************************************************/ /* Purpose: clement1_matrix() returns the CLEMENT1 matrix. Formula: if ( J = I + 1 ) A(I,J) = sqrt(I*(N-I)) else if ( I = J + 1 ) A(I,J) = sqrt(J*(N-J)) else A(I,J) = 0 Example: N = 5 . sqrt(4) . . . sqrt(4) . sqrt(6) . . . sqrt(6) . sqrt(6) . . . sqrt(6) . sqrt(4) . . . sqrt(4) . Properties: A is tridiagonal. A is banded, with bandwidth 3. Because A is tridiagonal, it has property A (bipartite). A is symmetric: A' = A. Because A is symmetric, it is normal. Because A is normal, it is diagonalizable. A is persymmetric: A(I,J) = A(N+1-J,N+1-I). The diagonal of A is zero. A is singular if N is odd. About 64 percent of the entries of the inverse of A are zero. The eigenvalues are plus and minus the numbers N-1, N-3, N-5, ..., (1 or 0). If N is even, det ( A ) = (-1)^(N/2) * (N-1) * (N+1)^(N/2) and if N is odd, det ( A ) = 0 Licensing: This code is distributed under the MIT license. Modified: 06 June 2024 Author: John Burkardt Reference: Paul Clement, A class of triple-diagonal matrices for test purposes, SIAM Review, Volume 1, 1959, pages 50-52. Input: int N, the order of the matrix. Output: double CLEMENT1[N*N], the matrix. */ { double *a; int i; int j; a = ( double * ) malloc ( n * n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { for ( j = 0; j < n; j++ ) { if ( j == i + 1 ) { a[i+j*n] = sqrt ( ( double ) ( ( i + 1 ) * ( n - i - 1 ) ) ); } else if ( i == j + 1 ) { a[i+j*n] = sqrt ( ( double ) ( ( j + 1 ) * ( n - j - 1 ) ) ); } else { a[i+j*n] = 0.0; } } } return a; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: timestamp() prints the current YMDHMS date as a time stamp. Example: 17 June 2014 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 01 May 2021 Author: John Burkardt */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); printf ( "%s\n", time_buffer ); return; # undef TIME_SIZE }