The linear system of equations Ax=0 or Ax = b , where A in an m×n matrix, x is an n×1 vector of unknowns and b is an m×1 constant vector can be solved by using the basic equation operations to reduce the equations to an upper triangular system or to an equation for the solution x. These equations can be represented by rows the matrix M=A or by the m×(n+1) augmented matrix matrix, M=(A,b). The three basic row operations are multiplying a row by a constant, swapping two rows or adding the multiple of one row to another.
The first program illustrates row operations. After entering a matrix, the user specifies the row operation and the program performs the desired operation.
/***********************************************************************************
Treibergs Mar. 14, 2006
Matrix Manipulator
Enter a matrix. Asks for desired row operation. Then does the operation.
Programming Notes
A type "matrix" is defined to be a structure that holds an array and two integers.
The main program asks for the desired operation, and calls the appropriate subroutine
to execute the operation. There are subroutines to enter the matrix, multiply one of
its rows by a constant, swap rows, add one row to another, and add a multiple of one
row to another.
Since C passes structure arguments "by value", assigning changes to the argument
within the subprogram do not affect the matrices in the calling program unless they
are returned. Therefore, we pass pointers to the structure to the subroutine.
(This is a bit tricky: see section 6.4 of Kernighan and Ritchie.) Thus in main,
"matrix a;" assigns memory to the variable a of type "matrix." The declaration
"matrix *a_ptr" assigns space for a single address to the pointer variable a_ptr,
and then the instruction "a_ptr = &a" puts the address of a in a_ptr.
Here is the advantage of passing pointers to the subprogram: pointers are passed
"by reference" so that the structure pointed to by the argument can be altered within
the subprogram. The function prototype is "void mulrom ( matrix *, int, double );"
which tells C that mulrom is a subprogram that takes a pointer to the type matrix,
an integer and a double argument and returns no value.
How to get the value of one of the variables in the structure matrix? For example,
the number of columns in the structure pointed at by a is "(*a_ptr).c". There is
another notation for this, namely, "a_ptr->c". Thus if we look inside the subprogram
void
mulrom ( matrix *a_ptr, int i, double c )
{
int j;
i=i-1;
for ( j=0; j < a_ptr->c; j++ )
(a_ptr->m)[i][j] = c * (a_ptr->m)[i][j];
return;
}
a_ptr is the local name of the pointer to data type matrix and "a_ptr->m[j][k]" or
"(a_ptr->m)[j][k]" is the value of the matrix at [j][k].
To contrast, look at the Gaussian Elimination program or Milicic's Row Operations
program to see how to pass arguments "by value." That way, you don't need to mess with
pointers and the mathematical logic is generally the same.
matmanip.c
*************************************************************************************/
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define MAXDIM 20
typedef struct
{
double m[MAXDIM][MAXDIM];
int r;
int c;
} matrix;
void
enterm ( matrix * );
void
entm( matrix * );
void
mulrom ( matrix *, int i, double c );
void
addrom ( matrix *, int i, int j );
void
addmurom ( matrix *, int i, double c, int j );
void
swaprom (matrix *, int i, int j );
void
printm ( matrix * );
int
main ( void )
{
double b;
int i=0, icount = 0, j, k, m, m1, m2;
matrix a;
matrix *a_ptr;
a_ptr = &a;
enterm ( a_ptr );
printm ( a_ptr );
while ( icount++ <= 1000 )
{
printf("What would you like to do? (1.) Multiply a row by a constant?\n"
" (2.) Divide a row by a constant?\n"
" (3.) Add one row to another?\n"
" (4.) Add a multiple of one row to another?\n"
" (5.) Swap two rows?\n"
" (6.) Enter a new matrix?\n"
" (7.) Quit?\n"
" Enter youir choice: (");
scanf ( "%d", &i );
switch ( i )
{
case 1:
printf ( "Row i is multiplied by c. Enter i c : ");
scanf ( "%d%lf", &j, &b );
if( ( 1 <=j ) && ( j <= a_ptr->r ) )
{
mulrom ( a_ptr, j, b );
printm ( a_ptr );
}
else
printf ( "Bad input: Row must be between %d and %d.\n", 1, a_ptr->r );
break;
case 2:
printf ( "Row i is divided by c. Enter i c : ");
scanf ( "%d%lf", &j, &b );
if( b == 0.0 )
printf ( "Bad input: Can't divide by zero.\n" );
else if( ( 1<=j ) && ( j <= a_ptr->r ) )
{
mulrom ( a_ptr, j, 1.0/b );
printm ( a_ptr );
}
else
printf ( "Bad input: Row must be between %d and %d.\n", 1, a_ptr->r );
break;
case 3:
printf ( "Row i is replaced by row i + row j. Enter i j : " );
scanf ( "%d%d", &j, &k );
if ( ( 1 <=j ) && ( j <= a_ptr->r ) && (1 <= k) && ( k <= a_ptr->r ) )
{
addrom ( a_ptr, j, k );
printm ( a_ptr );
}
else
printf ( "Bad input: Both rows must be between %d and %d.\n", 1, a_ptr->r );
break;
case 4:
printf ( "Row i is replaced by row i + c * row j. Enter i c j : ");
scanf ( "%d%lf%d", &j, &b, &k);
if ( ( 1 <=j ) && ( j <= a_ptr->r ) && ( 1 <= k ) && ( k <= a_ptr->c ) )
{
addmurom ( a_ptr , j, b, k );
printm ( a_ptr );
}
else
printf ( "Bad input: Both rows must be between %d and %d.\n", 1, a_ptr->r );
break;
case 5:
printf ( "Row i and row j are swapped. Enter i j : " );
scanf ( "%d%d", &j, &k );
if ( ( 1 <=j ) && ( j <= a_ptr->r ) && ( 1 <= k ) && ( k <= a_ptr->r ) )
{
swaprom ( a_ptr, j, k );
printm ( a_ptr );
}
else
printf ( "Bad input: Both rows must be between %d and %d.\n", 1, a_ptr->r );
break;
case 6:
enterm ( a_ptr );
printm ( a_ptr );
break;
case 7:
printf ( "Quit.\n" );
exit ( 0 );
break;
default:
break;
}
}
return EXIT_SUCCESS;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
enterm ( matrix *g_ptr )
{
int i=0, j, k, ir, ic;
while ( i++ <= 10 )
{
printf ( "\nEnter the number of rows and columns : " );
scanf ( "%d%d", &j, &k );
if ( (1 <=j) && (j <= MAXDIM) && (1 <= k) && (k <= MAXDIM) )
{
for ( ir=0; ir < j; ir++ )
{
for ( ic=0; ic < k; ic++ )
{
printf ( "Enter a(%d,%d) :", ir+1, ic+1 );
scanf ( "%lf", &g_ptr->m[ir][ic] );
}
}
g_ptr->r = j;
g_ptr->c = k;
return;
}
else
printf ( "Bad input: Both number of rows and number of columns "
"must be between %d and %d.\n", 1, MAXDIM );
}
exit ( 0 );
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
mulrom ( matrix *a_ptr, int i, double c )
{
int j;
i = i - 1;
for ( j=0; j < a_ptr->c; j++ )
(a_ptr->m)[i][j] = c * (a_ptr->m)[i][j];
return;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
addrom( matrix *a_ptr, int i, int j )
{
int k;
i = i - 1;
j = j - 1;
for ( k=0; k < a_ptr->c; k++ )
(a_ptr->m)[i][k] += (a_ptr->m)[j][k];
return;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
addmurom( matrix *a_ptr, int i, double c, int j )
{
int k;
i = i - 1;
j = j - 1;
for ( k=0; k < a_ptr->c; k++ )
(a_ptr->m)[i][k] += c * (a_ptr->m)[j][k];
return;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
swaprom (matrix *a_ptr, int i, int j)
{
int k;
double temp;
i = i - 1;
j = j - 1;
for ( k=0; k < a_ptr->c; k++ )
{
temp = a_ptr->m[i][k];
a_ptr->m[i][k] = a_ptr->m[j][k];
a_ptr->m[j][k] = temp;
}
return;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
printm ( matrix *a_ptr )
{
int i ,j;
printf ( "\n\n" );
for ( i=0; i < a_ptr->r; i++ )
{
for ( j=0; j < a_ptr->c; j++)
printf ( "\t%f", (a_ptr->m)[i][j] );
printf ( "\n\n" );
}
return;
}
/* Milicic
Elementary Row Transformations
*/
#include <stdio.h>
#include <math.h>
#define dim 20
typedef struct {
float coeff[dim][dim];
} matrix; /*define type matrix*/
matrix Matread (int n,int m)
{
int i, j;
matrix A;
printf("Enter matrix:\n");
for ( i=0 ; i < n ; i++ ) {
printf("\nEnter %d th row >> ",i+1);
for ( j=0 ; j < m ; j++ )
scanf("%f",&A.coeff[i][j]);
}
return(A);
}
/*Function Matdump dumps a matrix to the console */
void Matdump (matrix A, int n, int m)
{
int i, j;
for ( i=0 ; i < n ; i++ ) {
for ( j=0 ; j < m ; j++ )
printf("%f\t",A.coeff[i][j]);
printf("\n");
}
}
matrix Change (matrix A, int n, int m, int i, int j)
{
int k;
float x;
if (i != j) {
for ( k=0 ; k < m ; k++ ) {
x = A.coeff[i][k];
A.coeff[i][k] = A.coeff[j][k];
A.coeff[j][k] = x;
}
}
return(A);
}
/* Function AddRows adds ith row to jth row in the matrix A */
matrix AddRows (matrix A, int n, int m, int i, int j)
{
int k;
for ( k=0 ; k < m ; k++ )
A.coeff[j][k] = A.coeff[j][k] + A.coeff[i][k];
return(A);
}
/* Function MultRow multiplies elements of ith row in the matrix A by mult */
matrix MultRow (matrix A, int n, int m, int i, float mult)
{
int k;
for ( k=0 ; k < m ; k++ )
A.coeff[i][k] = mult * A.coeff[i][k];
return(A);
}
/* Function MultRow multiplies elements of ith row in the matrix A by mult
and adds them to the jth row */
matrix Rowmul (matrix A, int n, int m, int i, int j, float mult)
{
int k;
for ( k=0 ; k < m ; k++ )
A.coeff[j][k] = A.coeff[j][k] + A.coeff[i][k] * mult;
return(A);
}
int main() {
int n,m,i,j,menu,test;
float mult;
matrix A;
printf("Enter #rows, #columns >> ");
scanf("%d %d",&n, &m);
printf("\n");
A=Matread(n,m);
printf("\n");
Matdump(A,n,m);
test = 0 ;
do {
printf("\nDo you want to do a row transformation?\n");
printf("1. Switch two rows.\n");
printf("2. Add a row to another row.\n");
printf("3. Multiply a row by a number.\n");
printf("4. Multiply a row by a number and add it to another row.\n");
printf("0. Stop.\n");
printf("\nType your choice. >>> ");
scanf("%d",&menu);
switch(menu) {
case 1:
printf("\nWhich rows to switch? (i,j) >> ");
scanf("%d %d",&i, &j);
A=Change(A, n, m, i-1, j-1);
break;
case 2:
printf("\nAdd ith row to jth row (i,j) >> ");
scanf("%d %d",&i, &j);
A=AddRows(A, n, m, i-1, j-1);
break;
case 3:
printf("\nMultiply ith row with d (i,d) >> ");
scanf("%d %f",&i, &mult);
A=MultRow(A, n, m, i-1, mult);
break;
case 4:
printf("\nMultiply ith row with d and add it to jth row (i,j,d) >> ");
scanf("%d %d %f",&i, &j, &mult);
A=Rowmul(A, n, m, i-1, j-1, mult);
break;
default:
test = 1;
}
Matdump(A,n,m);
} while (test != 1);
printf("\n\nQuitting...\n");
}
The elmentary operations are done systematically to bring the matrix into row-echelon form and into reduced row-echelon form. From these reduced forms, the solutions of the linear system may be read off.
/***********************************************************************************
Treibergs Mar. 19, 2006
Gaussian Elimination
Enter a matrix. Finds the row-echelon and reduced-row-echelon form of the matrix.
Programming Notes
A type "matrix" is defined to be a structure that holds a double array and two
integers. The main program asks for matrix input. Then it does Gaussian Elimination
with partial pivoting to find the Echelon Form of the matrix. i.e., of the form
( 0 0 1 * * * * * * * * )
( 0 0 0 0 0 0 0 1 * * * )
( 0 0 0 0 0 0 0 0 1 * * )
( 0 0 0 0 0 0 0 0 0 0 0 )
where * denotes any number. This means that we go column by colums starting at the
leftmost (j=0) and from the first row lasti=0. From the jth row from (lasti=0)
to the last a.r, we find the row ipivot, that has the largest absolute value.
Then we swap the ipivot and ilast rows. Then we do elimination of all the rows i
from ilast to a.r below the pivotal row. This means divide the ilast row by
a[ilast][j] then subtract a[i][j] times ilast row from the ith row, leaving zeros
below the [ilast][j] element. Swapping in the largest value of the available
elements in the jth row to the ilast row is called partial pivoting. It
guarantees the best numerical performance, reducing the accuracy loss due to taking
the difference of near values which may be exaggerated when divideng by small values.
Finally, we call the routine "reducem" to eliminate the entries above the pivotal
elements. Thus the Reduced Row-Echelon form of the above matrix is
( 0 0 1 * * * * 0 0 * * )
( 0 0 0 0 0 0 0 1 0 * * )
( 0 0 0 0 0 0 0 0 1 * * )
( 0 0 0 0 0 0 0 0 0 0 0 )
Reducem assumes that the incoming matrix is in row-echelon form.
There are subroutines to enter the matrix, multiply one of its rows by a constant,
swap rows and add a multiple of one row to another. Since C passes structure
arguments "by value", the assigned changes to the argument within the subprogram
are made available to the the calling program via "return." The alternative is to
pass arguments "by reference.;" see my MatrixManip program.
Because of roundoff error, quantities that actually are zero may survive as
infinitessimal quantities. We therefore don't check that a quantity is zero but
rather whether the value is less than TOL which we've set to 0.5e-7. We assume
such small quantities are zero, but this may not in fact be true. It would mess
up calculations on really ill conditioned matrices.
gausselim.c
*************************************************************************************/
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define MAXDIM 20
# define TOL 0.5e-7
typedef struct
{
double m[MAXDIM][MAXDIM];
int r;
int c;
} matrix;
matrix
enterm ( void );
matrix
addmurom ( matrix, int i, double c, int j );
matrix
swaprom (matrix, int i, int j );
matrix
mulrom ( matrix a, int i, double c );
matrix
reducem ( matrix a );
void
printm ( matrix );
int
main ( void )
{
double pivot;
int i, j, k, ipivot, lasti=0, lastj=0;
matrix a;
a = enterm ();
printf("\n\nInitial Matrix");
printm ( a );
for( j=lastj; j < a.c; j++ )
{
ipivot=-1;
pivot=0.0;
/* Partial pivoting step: find the largest in pivotal col below last pivotal row */
for ( i = lasti; i<a.r; i++ )
{
if( fabs ( a.m[i][j] ) > fabs ( pivot ) + TOL )
{
pivot = a.m[i][j];
ipivot = i;
}
}
/* Skip to next col if all zeros. Else swap in pivoti row and eliminate other rows */
if ( ipivot != -1 )
{
if (ipivot != lasti)
a = swaprom ( a, ipivot, lasti);
if ( pivot != 1.0 )
a = mulrom ( a, lasti, 1.0/pivot );
for ( i = lasti + 1; i < a.r; i++ )
{
if ( fabs ( a.m[i][j] ) < TOL )
a.m[i][j] = 0.0;
else
a = addmurom ( a, i, -a.m[i][j], lasti );
}
lasti++;
}
}
printf("\nRow-Echelon Form of the Matrix");
printm ( a );
a = reducem ( a );
printf("\nReduced Row-Echelon Form of the Matrix");
printm ( a );
return EXIT_SUCCESS;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
matrix
enterm ( void )
{
int i=0, j, k, ir, ic;
matrix g;
while ( i++ <= 10 )
{
printf ( "\nEnter the number of rows and columns : " );
scanf ( "%d%d", &j, &k );
if ( (1<=j) && (j <= MAXDIM) && (1<= k) && (k <= MAXDIM) )
{
for ( ir=0; ir< j; ir++ )
{
for ( ic = 0; ic < k; ic++ )
{
printf ( "Enter a(%d,%d) :", ir+1, ic+1 );
scanf ( "%lf", &g.m[ir][ic] );
}
}
g.r = j;
g.c = k;
return g;
}
else
printf ( "Bad input: Both rows must be between %d and %d.\n", 1, MAXDIM );
}
exit ( 0 );
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
matrix
mulrom ( matrix a, int i, double c )
{
int j;
for ( j=0; j < a.c; j++ )
a.m[i][j] = c * a.m[i][j];
return a;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
matrix
addmurom( matrix a, int i, double c, int j )
{
int k;
for ( k = 0; k < a.c; k++ )
a.m[i][k] += c * a.m[j][k];
return a;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
matrix
swaprom (matrix a, int i, int j)
{
int k;
double temp;
for ( k=0; k < a.c; k++ )
{
temp = a.m[i][k];
a.m[i][k] = a.m[j][k];
a.m[j][k] = temp;
}
return a;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
printm ( matrix a )
{
int i, j;
printf ( "\n\n" );
for ( i = 0; i < a.r; i++ )
{
for ( j = 0; j < a.c; j++ )
printf ( "\t%f", a.m[i][j] );
printf ( "\n\n" );
}
return;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
matrix
reducem ( matrix a )
{
int i, j, ilast = 1, jlast = 1;
if ( a.r > 1)
{
while ( (ilast < a.r) && (jlast < a.c) )
{
if ( fabs ( a.m[ilast][jlast] ) < TOL )
jlast++;
else
{
for ( i=0; i<ilast; i++)
{
if ( fabs(a.m[i][jlast]) > TOL )
{
for ( j = jlast + 1; j<a.c; j++ )
a.m[i][j] -= a.m[i][jlast] * a.m[ilast][j];
}
a.m[i][jlast] = 0.0;
}
ilast++;
jlast++;
}
}
}
return a;
}
/* Milicic
Gauss Elimination
*/
#include <stdio.h>
#include <math.h>
#define true 1
#define false 0
#define dim 20
#define epsilon 0.000001
typedef struct {
float coeff[dim][dim];
} matrix; /*define type matrix*/
matrix Matread (int n,int m)
{
int i, j;
matrix A;
printf("Enter matrix:\n");
for ( i=0 ; i< n ; i++ ) {
printf("\nEnter %d th row >> ",i+1);
for ( j=0 ; j < m ; j++ )
scanf("%f",&A.coeff[i][j]);
}
return(A);
}
/*Function Matdump dumps a matrix to the console */
void Matdump (matrix A, int n, int m)
{
int i, j;
for ( i=0 ; i < n ; i++ ) {
for ( j=0 ; j < m ; j++ )
printf("%f\t",A.coeff[i][j]);
printf("\n");
}
}
matrix Change (matrix A, int n, int m, int i, int j)
{
int k;
float x;
if (i != j) {
for ( k=0 ; k < m ; k++ ) {
x = A.coeff[i][k];
A.coeff[i][k] = A.coeff[j][k];
A.coeff[j][k] = x;
}
}
return(A);
}
/* Function MultRow multiplies elements of ith row in the matrix A by mult */
matrix MultRow (matrix A, int n, int m, int i, float mult)
{
int k;
for ( k=0 ; k < m ; k++ )
A.coeff[i][k] = mult * A.coeff[i][k];
return(A);
}
/* Function Rowmul multiplies elements of ith row in the matrix A by mult
and adds them to the jth row */
matrix Rowmul (matrix A, int n, int m, int i, int j, float mult)
{
int k;
for ( k=0 ; k < m ; k++ )
A.coeff[j][k] = A.coeff[j][k] + A.coeff[i][k] * mult;
return(A);
}
/* Function Find finds the element of maximal absolute value in
the jth column below the element A (i,j) */
int Find(matrix A, int n, int m, int i, int j)
{
int k,p;
k = i;
for ( p = i+1 ; p < n ; p++ ) {
if (fabs(A.coeff[k][j]) < fabs(A.coeff[p][j]))
k = p;
}
return(k);
}
matrix Gauss (matrix A, int n, int m)
{
int i, j, counter;
int zero;
counter = 0;
for ( j=0 ; j < m ; j++ ) {
zero = true;
for ( i = j - counter ; i < n ; i++ )
if (fabs(A.coeff[i][j]) > epsilon)
zero = false;
if (zero == true)
counter = counter + 1;
else {
i=Find(A, n, m, j - counter,j);
A=Change(A, n, m, i, j - counter);
A=MultRow(A, n, m, j - counter,1 / A.coeff[j - counter][j]);
for ( i = j - counter + 1 ; i < n ; i++ )
if (fabs(A.coeff[i][j]) > epsilon)
A=Rowmul(A, n, m, j - counter, i, -A.coeff[i][j]);
}
}
return(A);
}
/* This procedure reduces a matrix A in row-echelon for to reduced
row-echelon form */
matrix Reduce (matrix A, int n, int m)
{
int i, j, k, zero;
for ( i = n-1 ; i > 0 ; i-- ) {
zero = true;
for ( j = 0 ; j < n ; j++ )
if (fabs(A.coeff[i][j]) > epsilon)
zero = false;
if (zero == false) {
j = 0;
while (fabs(A.coeff[i][j]) < epsilon)
j++;
for ( k = 0 ; k < i ; k++ )
if (fabs(A.coeff[k][j]) > epsilon)
A=Rowmul(A, n, m, i, k, -A.coeff[k][j]);
}
}
return(A);
}
int main() {
int n,m;
matrix A;
printf("Enter #rows, #columns >> ");
scanf("%d %d",&n, &m);
printf("\n");
A=Matread(n,m);
printf("\n");
Matdump(A, n, m);
A=Gauss(A, n, m);
printf("\n");
Matdump(A, n, m);
A=Reduce(A, n, m);
printf("\n");
Matdump(A, n, m);
}
/******************************************************
Treibergs 5-5-6
Multiplying matrices
Define the matrix structure.
Enter two matrices.
Check if compatible dimensions.
If so, multiply and print the result.
today.c
*******************************************************/
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define MAXDIM 20
typedef struct
{
double m[MAXDIM][MAXDIM];
int r;
int c;
}
matrix;
matrix
enterm ( void );
matrix
mulm ( matrix, matrix );
void
printm ( matrix );
int
main ( void )
{
matrix a, b, c;
printf("First matrix:\n");
a=enterm();
printf("\n\nSecond matrix:\n");
b=enterm();
c = mulm(a,b);
printf("\n\nFirst matrix:");
printm(a);
printf("\nSecond matrix:");
printm(b);
printf("\nProduct of the two matrices:");
printm(c);
return EXIT_SUCCESS;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
matrix
enterm ( void )
{
int i=0, j, k, ir, ic;
matrix g;
while ( i++ <= 10 )
{
printf ( "\nEnter the number of rows and columns : " );
scanf ( "%d%d", &j, &k );
if ( (1<=j) && (j <= MAXDIM) && (1<= k) && (k <= MAXDIM) )
{
for ( ir=0; ir< j; ir++ )
{
printf ( "Enter row a(%d,*) :", ir+1 );
for ( ic = 0; ic < k; ic++ )
scanf ( "%lf", &g.m[ir][ic] );
}
g.r = j;
g.c = k;
return g;
}
else
printf ( "Bad input: Both number of rows and number of columns must be between %d and %d.\n", 1, MAXDIM );
}
exit ( 0 );
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
matrix
mulm ( matrix a, matrix b )
{
double s;
int i,j,k;
matrix c;
if( a.c != b.r )
{
printf("Error--incompatible dimensions: the number of columns of the first matrix must\n"
" equal the number of rows of the second matrix.\n");
exit(0);
}
else
{
c.r=a.r;
c.c=b.c;
for( i=0; i<a.r; i++)
for(j=0; j<b.c; j++)
{
s=0.0;
for(k=0; k<b.r; k++)
s += a.m[i][k]*b.m[k][j];
c.m[i][j]=s;
}
}
return c;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
void
printm ( matrix a )
{
int i, j;
printf ( "\n\n" );
for ( i = 0; i < a.r; i++ )
{
for ( j = 0; j < a.c; j++ )
printf ( "\t%f", a.m[i][j] );
printf ( "\n\n" );
}
return;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */