matrix equation in C

[Tiza Naziri]

Hi,

How to represent this matrix equation in C code:

/ a0 \ / 0 1 0 1 0 0 1 0 \/ b0 \
| a1 | | 0 0 1 0 1 0 0 1 || b1 |
| a2 | | 1 0 0 1 0 1 0 0 || b2 |
| a3 |=| 0 1 0 0 1 0 1 0 || b3 |
| a4 | | 0 0 1 0 0 1 0 1 || b4 |
| a5 | | 1 0 0 1 0 0 1 0 || b5 |
| a6 | | 0 1 0 0 1 0 0 1 || b6 |
\ a7 / \ 1 0 1 0 0 1 0 0 /\ b7 /

Regards,

-- Tiza Naziri --

 

[Walter Roberson]

How to represent this matrix equation in C code:

[ 8 x 8 matrix multiplication shown ]

Depends -- do you need to do the calculation symbolically
or numerically? Are b0 thru b7 integers or floating point?
Are the matrix coefficients always going to be 0's and 1's with
no other possible values? Are you always going to be doing 8x8
or will the size be variable? Is the maximum size and density
enough that sparse representations seem like a good idea?
Are you solving for a[*] (matrix multiplication) or b[*] (matrix
division)? Is stability of the matrix inverse a significant concern?
Is there a good reason to impliment the manipulation yourself instead
of using one of the well-established libraries such as BLAS
or SCSL?

One hint: If you are doing the math yourself, then watch out for
the memory layout that C uses for multidimensional arrays, which
is different than the memory layout traditionally used by Fortran.
Pay attention to the memory layout when you are chosing which
index to loop over most quickly.

And if you get into very big matrices that might not fit in primary
cache, and you are implimenting by yourself, read up on matrix blocking
techniques.

Another hint: especially when your matrices have 2^N elements, watch
out for cache line aliasing -- if the array positions in memory are
seperated by a multiple of the cache line size, then something
seemingly simple such as

X[J] = Y[J];

might require a cache load to fetch the value from Y, and then might
require that the -same- cache line be dumped and loaded with the
corresponding position in X in order to do the store; whereas if
the starting positions of the arrays are carefully positioned to be
on different cache lines, then a cache-line's worth of data from Y
could be transfered to X before having to load cache again.
[Having a large primary cache doesn't help if the program's memory
access keeps thrashing over cache-line access...]


Anything to do with caches is outside the scope of the C89 standard.
Generally speaking pre-written libraries such as BLAS have been written
to take such matters into account, and so can sometimes end up being orders
of magnitude faster than naive code. It doesn't usually pay to
redevelop the mistakes of the past: it's usually best to find a
reputable matrix library that already does what you need.

A reference book for such matters: "Numerical Algorithms in C".

 

[Tiza Naziri]

Do you need to do the calculation symbolically
or numerically?
--> numerically..

Are b0 thru b7 integers or floating point?
--> b0 - b7 are integers of 0's and 1's..

Are the matrix coefficients always going to be 0's and 1's with
no other possible values?
--> yes (relates to previous question)..

Are you always going to be doing 8x8
or will the size be variable?
--> only 8x8 matrix

Are you solving for a[*] (matrix multiplication) or b[*] (matrix
division)?
--> a[*]

Any suggestion on this matter?

Regards.

-- Tiza Naziri --

 

[Rouben Rostamian]

How to represent this matrix equation in C code:

/ a0 \ / 0 1 0 1 0 0 1 0 \/ b0 \
| a1 | | 0 0 1 0 1 0 0 1 || b1 |
| a2 | | 1 0 0 1 0 1 0 0 || b2 |
| a3 |=| 0 1 0 0 1 0 1 0 || b3 |
| a4 | | 0 0 1 0 0 1 0 1 || b4 |
| a5 | | 1 0 0 1 0 0 1 0 || b5 |
| a6 | | 0 1 0 0 1 0 0 1 || b6 |
\ a7 / \ 1 0 1 0 0 1 0 0 /\ b7 /

How about this:

a0 = b1 + b3 + b6;
a1 = b2 + b4 + b7;
a2 = b0 + b3 + b5;
...

If this is not what you meant, then you will have
to be more explicit about your requirements.

 

[Walter Roberson]

--> b0 - b7 are integers of 0's and 1's..

Are the matrix coefficients always going to be 0's and 1's with

--> only 8x8 matrix

Are you solving for a[*] (matrix multiplication)

Any suggestion on this matter?

In the one special case of coefficients which are 0's and 1's,
multiplication is equivilent to bitwise AND, and you can carry
out the multiplications on groups at a time.

The below assumes that unsigned char is 8 bits or wider (which is
promised by the standard). If you needed to do larger matrices, you
could use wider datatypes than unsigned char.

#define MATSIZE 8

int sumbits( unsigned char v ) {
unsigned char t = v & 1;
unsigned char w = v >> 1;
int j;

for (j = 0; j < MATSIZE-1; j++ ) {
t += w & 1;
w >>= 1;
}

return t;
}

int main(void) {
unsigned char a[MATSIZE], M[MATSIZE];
unsigned char b;

init_M( M );
init_b( &b );

for (i = 0; i < MATSIZE; i++)
a = sumbits( M & b );

print_vec( a );
}

 

[Tiza Naziri]

Yes, but I want programming suggestions as written by walter.. Anyway,
thanks walter!