mat_mul_rect

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>

float *imsl_f_mat_mul_rect_band (char *string,

IMSL_A_MATRIX, int nrowa, int ncola, int nlca, int nuca, float a[],

IMSL_A_SYMMETRIC_STORAGE,

IMSL_B_MATRIX, int nrowb, int ncolb, int nlcb, int nucb, float b[],

IMSL_B_SYMMETRIC_STORAGE,

IMSL_X_VECTOR, int nx, float x[],

IMSL_RETURN_MATRIX_CODIAGONALS, int *nlc_result, int *nuc_result,

IMSL_RETURN_USER_VECTOR, float vector_user[],

Optional Arguments

IMSL_A_MATRIX, int nrowa, int ncola, int nlca, int nuca, float a[] (Input)
An array of size (nlca + 1 + nuca) × ncola containing the nrowa × ncola band matrix A with nlca lower and nuca upper codiagonals in band storage format.

IMSL_A_SYMMETRIC_STORAGE,

Indicates that matrix A is stored in band symmetric storage format, i.e., a is an array of size (nuca + 1) × ncola containing the diagonal and the nuca upper codiagonals of A.

IMSL_B_MATRIX, int nrowb, int ncolb, int nlcb, int nucb, float b[] (Input)
An array of size (nlcb + 1 + nucb) × ncolb containing the nrowb × ncolb band matrix B with nlcb lower and nucb upper codiagonals in band storage format.

IMSL_B_SYMMETRIC_STORAGE,

Indicates that matrix B is stored in band symmetric storage format, i.e., b is an array of size (nucb + 1) × ncolb containing the diagonal and the nucb upper codiagonals of B.

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 is a vector, return the answer in the user supplied array vector_user.

Description

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 M stands for A or B then the feasible matrix-vector products are "trans(x)*M", "trans(x)*trans(M)", "M*x" and "trans(M)*x". 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