mat_mul_rect_coordinate
Computes the transpose of a matrix, a matrix-vector product, or a matrix-matrix product for all matrices stored in sparse coordinate form.
Synopsis
#include <imsl.h>
void *imsl_f_mat_mul_rect_coordinate (char *string, ..., 0)
The equivalent double function is imsl_d_mat_mul_rect_coordinate.
Required Arguments
char *string (Input)
String indicating matrix multiplication to be performed.
Return Value
The returned value is the result of the multiplication. If the result is a vector, the return type is pointer to float. If the result of the multiplication is a sparse matrix, the return type is pointer to Imsl_f_sparse_elem. To release this space, use imsl_free.
Synopsis with Optional Arguments
#include <imsl.h>
void *imsl_f_mat_mul_rect_coordinate (char *string,
IMSL_A_MATRIX, int nrowa, int ncola, int nza, Imsl_f_sparse_elem *a,
IMSL_B_MATRIX, int nrowb, int ncolb, int nzb, Imsl_f_sparse_elem *b,
IMSL_X_VECTOR, int nx, float *x,
IMSL_RETURN_MATRIX_SIZE, int *size,
IMSL_RETURN_USER_VECTOR, float vector_user[],
0)
Optional Arguments
IMSL_A_MATRIX, int nrowa, int ncola, int nza, Imsl_f_sparse_elem *a (Input)
The sparse matrix
with nza nonzero elements.
IMSL_B_MATRIX, int nrowb, int ncolb, int nzb, Imsl_f_sparse_elem *b (Input)
The sparse matrix
with nzb nonzero elements.
IMSL_X_VECTOR, int nx, float *x, (Input)
The vector x of length nx.
IMSL_RETURN_MATRIX_SIZE, int *size, (Output)
If the function imsl_f_mat_mul_rect_coordinate returns a vector of type Imsl_f_sparse_elem, use this option to retrieve the length of the return vector, i.e. the number of nonzero elements in the sparse matrix generated by the requested computations.
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. It’s size depends on the computation.
Description
The function imsl_f_mat_mul_rect_coordinate computes a matrix-matrix product or a matrix-vector product, where the matrices are specified in coordinate representation. 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. Some multiplications, such as “A*trans(A)” or “trans(x)*B”, will produce a sparse matrix in coordinate format as a result. Other products such as “B*x” will produce a pointer to a floating type, containing the resulting vector.
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
In this example, a sparse matrix in coordinate form is multipled by a vector.
#include <imsl.h>
int main()
{
Imsl_f_sparse_elem a[] = {0, 0, 10.0,
1, 1, 10.0,
1, 2, -3.0,
1, 3, -1.0,
2, 2, 15.0,
3, 0, -2.0,
3, 3, 10.0,
3, 4, -1.0,
4, 0, -1.0,
4, 3, -5.0,
4, 4, 1.0,
4, 5, -3.0,
5, 0, -1.0,
5, 1, -2.0,
5, 5, 6.0};
float b[] = {10.0, 7.0, 45.0, 33.0, -34.0, 31.0};
int n = 6;
int nz = 15;
float *x;
/* Set x = A*b */
x = imsl_f_mat_mul_rect_coordinate ("A*x",
IMSL_A_MATRIX, n, n, nz, a,
IMSL_X_VECTOR, n, b,
0);
imsl_f_write_matrix ("Product Ab", 1, n, x, 0);
}
Output
Product Ab
1 2 3 4 5 6
100 -98 675 344 -302 162
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. 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 <math.h>
int main()
{
int i;
int n;
int c;
int nz;
int index;
Imsl_f_sparse_elem *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_coordinate (n, c, &nz, 0);
/* Convert to dense format */
for (i=0; i<nz; i++)
dense_a[a[i].col + n*a[i].row] = a[i].val;
/* 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_coordinate ("A*x",
IMSL_A_MATRIX, n, n, nz, 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.368035e-04