Can't get Jacobi algorithm to work in Objective-C

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 Gr,c and Gc,r. From the algorithm description, you want Gr,c=U1,2 (or u1.y) and Gc,r=U2,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 GT * A * G to 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.
Category:objective c Views:0 Time:2010-05-02

Related post

  • Algorithm to determine an objects neighbor 2009-02-05

    I vaguely remember reading about a programming exercise where objects are drawn on the screen. If an object has less than 2 neighbors, it dies because it is lonely, if it has more than 3 it dies because it is crowded. If the amount of neighbors is 2

  • JACOB doesn't release the objects properly 2009-06-11

    I have an eclipse plugin, which connects to a COM component using Jacob. But after I close the plugin entirely, the .exe file stays hanging in Windows processes. I use ComThread.InitMTA(true) for initialization and make sure that SafeRelease() is cal

  • implementing jacobi algorithm to implement laplace equation 2010-06-03

    The Algorithm traverses a 2D NxN array, making every element the average of its 4 surrounding neighbors (left, right, top, down). The NxN array has initially all zeros and is surrounded by a margin with all 1’s as shown in the example below. The 1’s

  • Algorithm for best positioning objects on a Visio model 2010-05-07

    I am trying to map all network devices and create a visio file with the resulting network topology. I was wondering if there are any algorithm for best positioning the nodes on the graph, considering its connections. Connections are bidirectional, li

  • Algorithm for finding nearest object on 2D grid 2010-07-25

    Say you have a 2D grid with each spot on the grid having x number of objects (with x >=0). I am having trouble thinking of a clean algorithm so that when a user specifies a coordinate, the algorithm finds the closest coordinate (including the one

  • Algorithm to Encrypt in Objective-C and Decrypt in Java 2010-11-09

    Want to know algorithm (if exists) that will match the results in both Objective-C and Java. Objective: Ojective-C will send the encrypted data and need to validate with the details that are encrypted in Java. So need the common algorithm where encry

  • C# Pathing algorithm for moving a object from point(X, Y) to point(X, Y) 2011-03-21

    Given a object which may move forward, backward, left and right at a given X,Y point. How to efficiently direct the object to a X,Y point using the given movement mechanics in the most efficient and human natural way. The object is available for move

  • quick pathfinding algorithm for lot of objects 2011-07-28

    i have lots of moving objects(maximum of 1000 objects) that each one needs to do pathfinding to maybe maximum 100 fixed locations(one object to one location at a time). What is the quickest pathfinding algorithm that best suit this? it doesn't have t

  • Algorithm to remove all objects of a tree from a list 2011-11-04

    I have a problem where I need to remove all objects of a tree from a list. I have a List<String> Tags which contains the tags in my entire system that match a certain criterion (generally starts with some search string). I also have a root Devi

  • Algorithm to generate combinatorial objects 2012-02-21

    I'm working on a timetable generation project and I need to pre-process the data before putting it into an LP model. I need to generate combinatorial objects to use in optimisation. The problem is very similar to the wood cutting problem. Say my I ha

  • Algorithm for sorting Java objects into buckets, then sorting within buckets 2012-04-30

    I have the need to take an ArrayList<Conference> conferences, where Conference contains a public Date beginDate parameter, and sort as follows: first, separate conferences into buckets which represent unique months for beginDate, and then sort

  • Algorithms for Tracking moving objects with a moving camera 2014-04-09

    I'm trying to develop an algorithm for real time tracking moving objects with a single moving camera setup as a project, in OpenCV (C++). My basic objectives are Detect motion in an (initially) static frame Track that moving object (camera to follow

  • Algorithm for randomly selecting object 2011-12-29

    I want to implement a simulation: there are 1000 objects; during a period of time 1800 seconds, each object is randomly selected (or whatever action); the number of selected objects along time follows a rough distribution: 30% will be selected within

  • adding force directed algorithm to Raphael SVG objects 2012-04-04

    I am building a user interface using Raphael JS. currently I have a .js script that draws out everything using Raphael JS 2.1 exactly as needed. However, because the drawing is driven by dynamic data it is highly likely that objects will overlap. Add

  • Red eye reduction algorithm 2008-09-25

    I need to implement red eye reduction for an application I am working on. Googling mostly provides links to commercial end-user products. Do you know a good red eye reduction algorithm, which could be used in a GPL application? --------------Solution

  • Approaching STL algorithms, lambda, local classes and other approaches 2008-11-02

    One of the things that seems to be necessary with use of STL is a way to specify local functions. Many of the functions that I would normally provide cannot be created using STL function object creation tools ( eg bind ), I have to hand roll my funct

  • Algorithm Speed Order of 2009-02-10

    Sometimes I get totally fooled trying to estimate an algorithm's speed with the O(x) notation, I mean, I can really point out when the order is O(n) or O(mxn), but for those that are O(lg(n)) or O(C(power n)) I think that I am missing something there

  • Algorithms for positioning objects in a canvas 2009-06-27

    Where can i find some algorithm to position some objects in canvas in a clever way? I'm using javascript (with Raphael svg library), but examples with other languages (or pseudo-language) are welcome. Geometry isn't my strong point =) For example hav

  • STL algorithms and const_iterators 2009-10-27

    Today I wrote a small predicate to find matching symbols in a container. But I'm faced to a problem: I want to use this predicate in a std::find_if call inside a const-method of a class, searching in a container that is a member of this class. But I

Copyright (C) dskims.com, All Rights Reserved.

processed in 0.139 (s). 11 q(s)