fork download
copy 
  1. std::vector<double> SquareMatrix::xVec(std::vector<double> bVec){
  2.  
  3. 	SquareMatrix temp = *this;
  4. 	std::vector<double> xVector(temp.dimension, 0);
  5.  
  6. 	//A^-1*b
  7. 	if (getDeterminant() != 0){
  8. 		SquareMatrix temp = *this;
  9. 		return temp.inverse()*bVec;
  10. 	}
  11. 	else if (temp.dimension != bVec.size()){
  12. 		std::cout << "Error: incompatible dimensions" << std::endl;
  13. 		return xVector;
  14. 	}
  15. 	else if (bVec.size() == 1 && bVec[0] == 0){
  16. 		xVector.resize(1);
  17. 		if (temp(1, 1) == 0){
  18. 			std::cout << "Any solution works." << std::endl;
  19. 		}
  20. 		xVector[0] = 0;
  21. 		return xVector;
  22. 	}
  23.  
  24. 	//row reduces the matrix temp.
  25. 	for (int i = 0; i < temp.dimension - 1; ++i){
  26.  
  27. 		std::cout << "First, at i: " << i << std::endl;
  28. 		temp.show();
  29.  
  30. 		//swaps the rows if pivot point is zero
  31. 		if (temp.matrix[i][i] == 0){
  32. 			for (int k = i + 1; k < temp.dimension; ++k){
  33. 				if (temp.matrix[k][i] != 0){
  34. 					temp.rowSwap(k + 1, i + 1);
  35. 					std::swap(bVec[k], bVec[i]);
  36. 					break;
  37. 				}
  38. 			}
  39. 		}
  40.  
  41. 		std::cout << "Next: " << std::endl;
  42. 		temp.show();
  43.  
  44. 		//multiplies entire row by reciprocal of pivot.
  45. 		for (int j = i; j < temp.dimension; ++j){
  46. 			if (temp.matrix[j][i] != 0){
  47. 				//bVec[j] reassignment MUST go before
  48. 				//temp reassignment.
  49. 				bVec[j] /= temp.matrix[j][i];
  50. 				temp.rowMul(j + 1, 1 / temp.matrix[j][i]);
  51. 			}
  52. 		}
  53.  
  54. 		std::cout << "Next: " << std::endl;
  55. 		temp.show();
  56.  
  57. 		//subtracts rows from each other
  58. 		for (int j = i + 1; j < temp.dimension; ++j){
  59. 			if (temp.matrix[j][i] != 0){
  60. 				temp.rowNeg(j + 1, i + 1);
  61. 				bVec[j] -= bVec[i];
  62. 			}
  63. 		}
  64.  
  65. 		std::cout << "Next: " << std::endl;
  66. 		temp.show();
  67.  
  68. 	}
  69.  
  70. 	if (temp.matrix[dimension - 1][dimension - 1] != 0){
  71. 		bVec[dimension - 1] /= temp.matrix[dimension - 1][dimension - 1];
  72. 		temp.matrix[dimension - 1][dimension - 1] 
  73. 			/= temp.matrix[dimension - 1][dimension - 1];
  74. 	}
  75.  
  76. 	//(inefficient) bubble-sort algorithm to 
  77. 	//place the pivots where they belong
  78. 	for (int i = 0; i < temp.dimension; ++i){
  79. 		for (int j = 0; j < dimension; j++){
  80. 			if (temp.matrix[i][j] != 1){
  81. 				continue;
  82. 			}
  83. 			else if (i != j){
  84. 				temp.rowSwap(i + 1, j + 1);
  85. 			}
  86. 		}
  87. 	}
  88.  
  89. 	//matrix A is row-reduced and b is consistent with these operations
  90. 	//now we are ready for back substitution
  91.  
  92. 	//the last index
  93. 	int start = temp.dimension - 1;
  94.  
  95. 	for (int k = 0; k < temp.dimension; ++k){
  96. 		double restOfRowValue = 0;
  97. 		int counter = 0;
  98. 		//detects if equation has a free variable.
  99. 		if (bVec[start - k] == 0){
  100. 			for (int i = 0; i < temp.dimension; ++i){
  101. 				if (temp.matrix[start - k][i] == 0){
  102. 					++counter;
  103. 				}
  104. 				else
  105. 					break;
  106. 			}
  107. 			if (counter == temp.dimension){
  108. 				std::cout << "x_" << start - k + 1 << " is free;";
  109. 				std::cout << " setting it equal to one..." << std::endl;
  110. 				xVector[start - k] = 1;
  111. 				continue;
  112. 			}
  113. 		}
  114.  
  115. 		//detects if equation is inconsistent
  116. 		else if (temp.matrix[start - k][start - k] == 0){
  117. 			std::cout << "Inconsistent system detected;";
  118. 			std::cout << " matrix has only trivial solution." << std::endl;
  119. 			std::vector<double> newVec(temp.dimension, 0);
  120. 			return newVec;
  121. 		}
  122.  
  123. 		//b - Ax to find the next x-value
  124. 		for (int j = start; start - k < j; --j){
  125. 			restOfRowValue += temp.matrix[start - k][j]*xVector[j];
  126. 		}
  127.  
  128. 		xVector[start - k] = bVec[start - k] - restOfRowValue;
  129.  
  130. 	}
  131.  
  132. 	return xVector;
  133.  
  134. }
Compilation error #stdin compilation error #stdout 0s 0KB
comments (?)
stdin
copy 
Standard input is empty
compilation info 
prog.cpp:1:6: error: 'vector' in namespace 'std' does not name a template type
 std::vector<double> SquareMatrix::xVec(std::vector<double> bVec){
      ^
stdout
copy 
Standard output is empty
https://ideone.com/oCefy9
language:
C++14 (gcc 8.3)
created:
10 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!