Can't get Jacobi algorithm to work in Objective-C
- by Chris Long
Hi,
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 in advance,
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++ ) {
value = 0.0;
for ( j = 0; j < 5; j++ ) {
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];
}
// ...
