For some reason, I can't get this program to work. I've had other CS majors look at it and they can't figure it out either.

This program performs the Jacobi algorithm (you can see step-by-step instructions and a MATLAB implementation here). BTW, it's different from the Wikipedia article of the same name.

Since `NSArray`

is one-dimensional, I added a method that makes it act like a two-dimensional C array. After running the Jacobi algorithm many times, the diagonal entries in the `NSArray`

(`i[0][0]`

, `i[1][1]`

, etc.) are supposed to get bigger and the others approach 0. For some reason though, they all increase exponentially. For instance, `i[2][4]`

should equal `0.0000009`

, not `9999999`

, while `i[2][2]`

should be big.

Thanks,

Chris

##
NSArray+Matrix.m

`@implementation NSArray (Matrix) @dynamic offValue, transposed; - (double)offValue { double sum = 0.0; for ( MatrixItem *item in self ) if ( item.nonDiagonal ) sum += pow( item.value, 2.0 ); return sum; } - (NSMutableArray *)transposed { NSMutableArray *transpose = [[[NSMutableArray alloc] init] autorelease]; int i, j; for ( i = 0; i < 5; i++ ) { for ( j = 0; j < 5; j++ ) { [transpose addObject:[self objectAtRow:j andColumn:i]]; } } return transpose; } - (id)objectAtRow:(NSUInteger)row andColumn:(NSUInteger)column { NSUInteger index = 5 * row + column; return [self objectAtIndex:index]; } - (NSMutableArray *)multiplyWithMatrix:(NSArray *)array { NSMutableArray *result = [[NSMutableArray alloc] init]; int i = 0, j = 0, k = 0; double value; for ( i = 0; i < 5; i++ ) { for ( j = 0; j < 5; j++ ) { value = 0.0; // (JeremyP's answer) for ( k = 0; k < 5; k++ ) { MatrixItem *firstItem = [self objectAtRow:i andColumn:k]; MatrixItem *secondItem = [array objectAtRow:k andColumn:j]; value += firstItem.value * secondItem.value; } MatrixItem *item = [[MatrixItem alloc] initWithValue:value]; item.row = i; item.column = j; [result addObject:item]; } } return result; } @end `

##
Jacobi_AlgorithmAppDelegate.m

`// ... - (void)jacobiAlgorithmWithEntry:(MatrixItem *)entry { MatrixItem *b11 = [matrix objectAtRow:entry.row andColumn:entry.row]; MatrixItem *b22 = [matrix objectAtRow:entry.column andColumn:entry.column]; double muPlus = ( b22.value + b11.value ) / 2.0; muPlus += sqrt( pow((b22.value - b11.value), 2.0) + 4.0 * pow(entry.value, 2.0) ); Vector *u1 = [[[Vector alloc] initWithX:(-1.0 * entry.value) andY:(b11.value - muPlus)] autorelease]; [u1 normalize]; Vector *u2 = [[[Vector alloc] initWithX:-u1.y andY:u1.x] autorelease]; NSMutableArray *g = [[[NSMutableArray alloc] init] autorelease]; for ( int i = 0; i <= 24; i++ ) { MatrixItem *item = [[[MatrixItem alloc] init] autorelease]; if ( i == 6*entry.row ) item.value = u1.x; else if ( i == 6*entry.column ) item.value = u2.y; else if ( i == ( 5*entry.row + entry.column ) || i == ( 5*entry.column + entry.row ) ) item.value = u1.y; else if ( i % 6 == 0 ) item.value = 1.0; else item.value = 0.0; [g addObject:item]; } NSMutableArray *firstResult = [[g.transposed multiplyWithMatrix:matrix] autorelease]; matrix = [firstResult multiplyWithMatrix:g]; } // ... `

-------------Problems Reply------------

When you add the square root term to `muPlus`

, you don't divide by two. The calculation should be either:

`double muPlus = ( b22.value + b11.value ) / 2.0;`

muPlus += sqrt( pow((b22.value - b11.value), 2.0)

+ 4.0 * pow(entry.value, 2.0)

) / 2.0;

or:

`double muPlus = ( b22.value + b11.value );`

muPlus += sqrt( pow((b22.value - b11.value), 2.0)

+ 4.0 * pow(entry.value, 2.0) );

muPlus /= 2.0;

Also, you assign `u1.y`

to both G_{r,c} and G_{c,r}. From the algorithm description, you want G_{r,c}=U_{1,2} (or `u1.y`

) and G_{c,r}=U_{2,1} (or `u2.x`

). Note that you don't actually need `u2`

; you can substitute `-u1.y`

for `u2.x`

and `u1.x`

for `u2.y`

.

###
Off-Topic

According to the Fundamental Rule of Cocoa Memory Management, `-[NSArray multiplyWithMatrix:]`

should return an autoreleased array, since the multiplicand should relinquish ownership. Also, you should use accessors to assign `G`

to ^{T} * A * G`matrix`

rather than doing it directly so that it can be properly managed.

Since most of the tests in the loop to fill out `g`

will be false during each iteration, it's most likely more efficient to fill `g`

with some default values and then update `g`

. You could create a zero matrix, then set the diagonal to ones, then fill in the values from U, or you could create an identity matrix (leave the `i%6 == 0`

test in the loop) then fill in the values from U. Profile each of the three approaches.

Have you got any unit tests for your matrix category? I mean, are you certain that the multiplication algorithm works? I would say that initialising value to 0 happens in the wrong loop. I think you need to do it inside the j loop.

A couple of other observations:

- You don't need the
`@dynamic`

property declaration because you are defining the implementation of the properties yourself. - Consider creating your own Matrix class that wraps a normal C array of doubles. You might find the implementation a bit simpler.