mat_mul_rect_band
Computes the transpose of a matrix, a matrix-vector product, or a matrix-matrix product, all matrices stored in band form.
Synopsis
#include <imsl.h>
float *imsl_f_mat_mul_rect_band (char *string, ..., 0)
The equivalent double function is imsl_d_mat_mul_rect_band.
Required Arguments
char *string (Input)
String indicating matrix multiplication to be performed.
Return Value
The result of the multiplication is returned. To release this space, use imsl_free.
Synopsis with Optional Arguments
#include <imsl.h>
void *imsl_f_mat_mul_rect_band (char *string,
IMSL_A_MATRIX, int nrowa, int ncola, int nlca, int nuca, float *a,
IMSL_B_MATRIX, int nrowb, int ncolb, int nlcb, int nucb, float *b,
IMSL_X_VECTOR, int nx, float *x,
IMSL_RETURN_MATRIX_CODIAGONALS, int *nlc_result, int *nuc_result,
IMSL_RETURN_USER_VECTOR, float vector_user[],
0)
Optional Arguments
IMSL_A_MATRIX, int nrowa, int ncola, int nlca, int nuca, float *a (Input)
The sparse matrix
IMSL_B_MATRIX, int nrowb, int ncolb, int nlcb, int nucb, float *b (Input)
The sparse matrix
IMSL_X_VECTOR, int nx, float *x, (Input)
The vector x of length nx.
IMSL_RETURN_MATRIX_CODIAGONALS, int *nlc_result, int *nuc_result, (Output)
If the function imsl_f_mat_mul_rect_band returns data for a band matrix, use this option to retrieve the number of lower and upper codiagonals of the return matrix.
IMSL_RETURN_USER_VECTOR, float vector_user[], (Output)
If the result of the computation in a vector, return the answer in the user supplied sparse vector_user.
Description
The function imsl_f_mat_mul_rect_band computes a matrix-matrix product or a matrix-vector product, where the matrices are specified in band format. The operation performed is specified by string. For example, if “A*x” is given, Ax is computed. In string, the matrices A and B and the vector x can be used. Any of these names can be used with trans, indicating transpose. The vector x is treated as a dense n × 1 matrix. If string contains only one item, such as “x” or “trans(A)”, then a copy of the array, or its transpose is returned.
The matrices and/or vector referred to in string must be given as optional arguments. Therefore, if string is “A*x”, then IMSL_A_MATRIX and IMSL_X_VECTOR must be given.
Examples
Example 1
Consider the matrix
After storing A in band format, multiply A by x = (1, 2, 3, 4)T and print the result.
#include <imsl.h>
int main()
{
float a[] = {0.0, -1.0, -2.0, 2.0,
2.0, 1.0, -1.0, 1.0,
-3.0, 0.0, 2.0, 0.0};
float x[] = {1.0, 2.0, 3.0, 4.0};
int n = 4;
int nuca = 1;
int nlca = 1;
float *b;
/* Set b = A*x */
b = imsl_f_mat_mul_rect_band ("A*x",
IMSL_A_MATRIX, n, n, nlca, nuca, a,
IMSL_X_VECTOR, n, x,
0);
imsl_f_write_matrix ("Product, Ax", 1, n, b, 0);
}
Output
Product, Ax
1 2 3 4
0 -7 5 10
Example 2
This example uses the power method to determine the dominant eigenvector of E(100, 10). The same computation is performed by using imsl_f_eig_sym, described in the chapter Eigensystem Analysis. The iteration stops when the component-wise absolute difference between the dominant eigenvector found by imsl_f_eig_sym and the eigenvector at the current iteration is less than the square root of machine unit roundoff.
#include <imsl.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main()
{
int i;
int j;
int k;
int n;
int c;
int nz;
int index;
int start;
int stop;
float *a;
float *z;
float *q;
float *dense_a;
float *dense_evec;
float *dense_eval;
float norm;
float *evec;
float error;
float tolerance;
n = 100;
c = 10;
tolerance = sqrt(imsl_f_machine(4));
error = 1.0;
evec = (float*) malloc (n*sizeof(*evec));
z = (float*) malloc (n*sizeof(*z));
q = (float*) malloc (n*sizeof(*q));
dense_a = (float*) calloc (n*n, sizeof(*dense_a));
a = imsl_f_generate_test_band (n, c,
0);
/* Convert to dense format,
starting with upper triangle */
start = c;
for (i=0; i<c; i++, start--)
for (k=0, j=start; j<n; j++, k++)
dense_a[k*n + j] = a[i*n + j];
/* Convert diagonal */
for (j=0; j<n; j++)
dense_a[j*n + j] = a[c*n + j];
/* Convert lower triangle */
stop = n-1;
for (i=c+1; i<2*c+1; i++, stop--)
for (k=i-c, j=0; j<stop; j++, k++)
dense_a[k*n + j] = a[i*n + j];
/* Determine dominant eigenvector by a dense method */
dense_eval = imsl_f_eig_sym (n, dense_a,
IMSL_VECTORS, &dense_evec,
0);
for (i=0; i<n; i++)
evec[i] = dense_evec[n*i];
/* Normalize */
norm = imsl_f_vector_norm (n, evec,
0);
for (i=0; i<n; i++)
evec[i] /= norm;
for (i=0; i<n; i++)
q[i] = 1.0/sqrt((float) n);
/* Do power method */
while (error > tolerance) {
imsl_f_mat_mul_rect_band ("A*x",
IMSL_A_MATRIX, n, n, c, c, a,
IMSL_X_VECTOR, n, q,
IMSL_RETURN_USER_VECTOR, z,
0);
/* Normalize */
norm = imsl_f_vector_norm (n, z,
0);
for (i=0; i<n; i++)
q[i] = z[i]/norm;
/* Compute maximum absolute error between any
two elements */
error = imsl_f_vector_norm (n, q,
IMSL_SECOND_VECTOR, evec,
IMSL_INF_NORM, &index,
0);
}
printf ("Maximum absolute error = %e\n", error);
}
Output
Maximum absolute error = 3.367960e-04
