circuitjs1 icon indicating copy to clipboard operation
circuitjs1 copied to clipboard

How do I feel that the lu_factor method in the CirSim class is Doolittle decomposition?

Open 306Vulcan opened this issue 2 years ago • 0 comments

Hello everyone, can you help me take a look. I looked at the code and divided it by the pivot when calculating the lower triangular part, which seems to be Doolittle decomposition. Why Crout decomposition?

` // factors a matrix into upper and lower triangular matrices by // gaussian elimination. On entry, a[0..n-1][0..n-1] is the // matrix to be factored. ipvt[] returns an integer vector of pivot // indices, used in the lu_solve() routine. static boolean lu_factor(double a[][], int n, int ipvt[]) { int i, j, k;

	// check for a possible singular matrix by scanning for rows that
	// are all zeroes
	for (i = 0; i != n; i++) {
		boolean row_all_zeros = true;
		for (j = 0; j != n; j++) {
			if (a[i][j] != 0) {
				row_all_zeros = false;
				break;
			}
		}
		// if all zeros, it's a singular matrix
		if (row_all_zeros)
			return false;
	}

	// use Crout's method; loop through the columns
	for (j = 0; j != n; j++) {

		// calculate upper triangular elements for this column
		for (i = 0; i != j; i++) {
			double q = a[i][j];
			for (k = 0; k != i; k++)
				q -= a[i][k] * a[k][j];
			a[i][j] = q;
		}

		// calculate lower triangular elements for this column
		double largest = 0;
		int largestRow = -1;
		for (i = j; i != n; i++) {
			double q = a[i][j];
			for (k = 0; k != j; k++)
				q -= a[i][k] * a[k][j];
			a[i][j] = q;
			double x = Math.abs(q);
			if (x >= largest) {
				largest = x;
				largestRow = i;
			}
		}

		// pivoting
		if (j != largestRow) {
			double x;
			for (k = 0; k != n; k++) {
				x = a[largestRow][k];
				a[largestRow][k] = a[j][k];
				a[j][k] = x;
			}
		}

		// keep track of row interchanges
		ipvt[j] = largestRow;

		// avoid zeros
		if (a[j][j] == 0.0) {
			System.out.println("avoided zero");
			a[j][j] = 1e-18;
		}

		if (j != n - 1) {
			double mult = 1.0 / a[j][j];
			for (i = j + 1; i != n; i++)
				a[i][j] *= mult;
		}
	}
	return true;
}`

306Vulcan avatar Apr 12 '23 03:04 306Vulcan