PageRenderTime 63ms CodeModel.GetById 32ms RepoModel.GetById 1ms app.codeStats 0ms

/Utilities/vxl/v3p/netlib/sparse/spOutput.c

https://github.com/luisibanez/ITK
C | 378 lines | 205 code | 39 blank | 134 comment | 60 complexity | 037519c2d9e6c21064ceb917a7851173 MD5 | raw file
  1. /*
  2. * MATRIX OUTPUT MODULE
  3. *
  4. * Author: Advisor:
  5. * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
  6. * UC Berkeley
  7. */
  8. /*! \file
  9. *
  10. * This file contains the output-to-file and output-to-screen routines for
  11. * the matrix package.
  12. *
  13. * Objects that begin with the \a spc prefix are considered private
  14. * and should not be used.
  15. *
  16. * \author
  17. * Kenneth S. Kundert <kundert@users.sourceforge.net>
  18. */
  19. /* >>> User accessible functions contained in this file:
  20. * spPrint
  21. * spFileMatrix
  22. * spFileVector
  23. * spFileStats
  24. *
  25. * >>> Other functions contained in this file:
  26. */
  27. /*
  28. * Revision and copyright information.
  29. *
  30. * Copyright (c) 1985-2003
  31. * by Kenneth S. Kundert
  32. */
  33. #if 0
  34. static char copyright[] =
  35. "Sparse1.4: Copyright (c) 1985-2003 by Kenneth S. Kundert";
  36. #endif
  37. /*
  38. Removed File IO routines to get rid of fopen warnings - JLM
  39. */
  40. /*
  41. * IMPORTS
  42. *
  43. * >>> Import descriptions:
  44. * spConfig.h
  45. * Macros that customize the sparse matrix routines.
  46. * spMatrix.h
  47. * Macros and declarations to be imported by the user.
  48. * spDefs.h
  49. * Matrix type and macro definitions for the sparse matrix routines.
  50. */
  51. #define spINSIDE_SPARSE
  52. #include <stdio.h>
  53. #include "spConfig.h"
  54. #include "spMatrix.h"
  55. #include "spDefs.h"
  56. #if DOCUMENTATION
  57. /*!
  58. * Formats and send the matrix to standard output. Some elementary
  59. * statistics are also output. The matrix is output in a format that is
  60. * readable by people.
  61. *
  62. * \param eMatrix
  63. * Pointer to matrix.
  64. * \param PrintReordered
  65. * Indicates whether the matrix should be printed out in its original
  66. * form, as input by the user, or whether it should be printed in its
  67. * reordered form, as used by the matrix routines. A zero indicates that
  68. * the matrix should be printed as inputed, a one indicates that it
  69. * should be printed reordered.
  70. * \param Data
  71. * Boolean flag that when false indicates that output should be
  72. * compressed such that only the existence of an element should be
  73. * indicated rather than giving the actual value. Thus 11 times as
  74. * many can be printed on a row. A zero signifies that the matrix
  75. * should be printed compressed. A one indicates that the matrix
  76. * should be printed in all its glory.
  77. * \param Header
  78. * Flag indicating that extra information should be given, such as row
  79. * and column numbers.
  80. */
  81. /* >>> Local variables:
  82. * Col (int)
  83. * Column being printed.
  84. * ElementCount (int)
  85. * Variable used to count the number of nonzero elements in the matrix.
  86. * LargestElement (RealNumber)
  87. * The magnitude of the largest element in the matrix.
  88. * LargestDiag (RealNumber)
  89. * The magnitude of the largest diagonal in the matrix.
  90. * Magnitude (RealNumber)
  91. * The absolute value of the matrix element being printed.
  92. * PrintOrdToIntColMap (int [])
  93. * A translation array that maps the order that columns will be
  94. * printed in (if not PrintReordered) to the internal column numbers.
  95. * PrintOrdToIntRowMap (int [])
  96. * A translation array that maps the order that rows will be
  97. * printed in (if not PrintReordered) to the internal row numbers.
  98. * pElement (ElementPtr)
  99. * Pointer to the element in the matrix that is to be printed.
  100. * pImagElements (ElementPtr [ ])
  101. * Array of pointers to elements in the matrix. These pointers point
  102. * to the elements whose real values have just been printed. They are
  103. * used to quickly access those same elements so their imaginary values
  104. * can be printed.
  105. * Row (int)
  106. * Row being printed.
  107. * Size (int)
  108. * The size of the matrix.
  109. * SmallestDiag (RealNumber)
  110. * The magnitude of the smallest diagonal in the matrix.
  111. * SmallestElement (RealNumber)
  112. * The magnitude of the smallest element in the matrix excluding zero
  113. * elements.
  114. * StartCol (int)
  115. * The column number of the first column to be printed in the group of
  116. * columns currently being printed.
  117. * StopCol (int)
  118. * The column number of the last column to be printed in the group of
  119. * columns currently being printed.
  120. * Top (int)
  121. * The largest expected external row or column number.
  122. */
  123. void
  124. spPrint(
  125. spMatrix eMatrix,
  126. int PrintReordered,
  127. int Data,
  128. int Header
  129. )
  130. {
  131. MatrixPtr Matrix = (MatrixPtr)eMatrix;
  132. register int J = 0;
  133. int I, Row, Col, Size, Top, StartCol = 1, StopCol, Columns, ElementCount = 0;
  134. double Magnitude, SmallestDiag, SmallestElement;
  135. double LargestElement = 0.0, LargestDiag = 0.0;
  136. ElementPtr pElement, pImagElements[PRINTER_WIDTH/10+1];
  137. int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
  138. /* Begin `spPrint'. */
  139. ASSERT_IS_SPARSE( Matrix );
  140. Size = Matrix->Size;
  141. /* Create a packed external to internal row and column translation array. */
  142. # if TRANSLATE
  143. Top = Matrix->AllocatedExtSize;
  144. #else
  145. Top = Matrix->AllocatedSize;
  146. #endif
  147. CALLOC( PrintOrdToIntRowMap, int, Top + 1 );
  148. CALLOC( PrintOrdToIntColMap, int, Top + 1 );
  149. if ( PrintOrdToIntRowMap == NULL OR PrintOrdToIntColMap == NULL)
  150. { Matrix->Error = spNO_MEMORY;
  151. return;
  152. }
  153. for (I = 1; I <= Size; I++)
  154. { PrintOrdToIntRowMap[ Matrix->IntToExtRowMap[I] ] = I;
  155. PrintOrdToIntColMap[ Matrix->IntToExtColMap[I] ] = I;
  156. }
  157. /* Pack the arrays. */
  158. for (J = 1, I = 1; I <= Top; I++)
  159. { if (PrintOrdToIntRowMap[I] != 0)
  160. PrintOrdToIntRowMap[ J++ ] = PrintOrdToIntRowMap[ I ];
  161. }
  162. for (J = 1, I = 1; I <= Top; I++)
  163. { if (PrintOrdToIntColMap[I] != 0)
  164. PrintOrdToIntColMap[ J++ ] = PrintOrdToIntColMap[ I ];
  165. }
  166. /* Print header. */
  167. if (Header)
  168. { printf("MATRIX SUMMARY\n\n");
  169. printf("Size of matrix = %1d x %1d.\n", Size, Size);
  170. if ( Matrix->Reordered AND PrintReordered )
  171. printf("Matrix has been reordered.\n");
  172. putchar('\n');
  173. if ( Matrix->Factored )
  174. printf("Matrix after factorization:\n");
  175. else
  176. printf("Matrix before factorization:\n");
  177. SmallestElement = LARGEST_REAL;
  178. SmallestDiag = SmallestElement;
  179. }
  180. if (Size == 0) return;
  181. /* Determine how many columns to use. */
  182. Columns = PRINTER_WIDTH;
  183. if (Header) Columns -= 5;
  184. if (Data) Columns = (Columns+1) / 10;
  185. /*
  186. * Print matrix by printing groups of complete columns until all the columns
  187. * are printed.
  188. */
  189. J = 0;
  190. while ( J <= Size )
  191. /* Calculate index of last column to printed in this group. */
  192. { StopCol = StartCol + Columns - 1;
  193. if (StopCol > Size)
  194. StopCol = Size;
  195. /* Label the columns. */
  196. if (Header)
  197. { if (Data)
  198. { printf(" ");
  199. for (I = StartCol; I <= StopCol; I++)
  200. { if (PrintReordered)
  201. Col = I;
  202. else
  203. Col = PrintOrdToIntColMap[I];
  204. printf(" %9d", Matrix->IntToExtColMap[ Col ]);
  205. }
  206. printf("\n\n");
  207. }
  208. else
  209. { if (PrintReordered)
  210. printf("Columns %1d to %1d.\n",StartCol,StopCol);
  211. else
  212. { printf("Columns %1d to %1d.\n",
  213. Matrix->IntToExtColMap[ PrintOrdToIntColMap[StartCol] ],
  214. Matrix->IntToExtColMap[ PrintOrdToIntColMap[StopCol] ]);
  215. }
  216. }
  217. }
  218. /* Print every row ... */
  219. for (I = 1; I <= Size; I++)
  220. { if (PrintReordered)
  221. Row = I;
  222. else
  223. Row = PrintOrdToIntRowMap[I];
  224. if (Header)
  225. { if (PrintReordered AND NOT Data)
  226. printf("%4d", I);
  227. else
  228. printf("%4d", Matrix->IntToExtRowMap[ Row ]);
  229. if (NOT Data) putchar(' ');
  230. }
  231. /* ... in each column of the group. */
  232. for (J = StartCol; J <= StopCol; J++)
  233. { if (PrintReordered)
  234. Col = J;
  235. else
  236. Col = PrintOrdToIntColMap[J];
  237. pElement = Matrix->FirstInCol[Col];
  238. while(pElement != NULL AND pElement->Row != Row)
  239. pElement = pElement->NextInCol;
  240. if (Data)
  241. pImagElements[J - StartCol] = pElement;
  242. if (pElement != NULL)
  243. /* Case where element exists */
  244. { if (Data)
  245. printf(" %9.3g", (double)pElement->Real);
  246. else
  247. putchar('x');
  248. /* Update status variables */
  249. if ( (Magnitude = ELEMENT_MAG(pElement)) > LargestElement )
  250. LargestElement = Magnitude;
  251. if ((Magnitude < SmallestElement) AND (Magnitude != 0.0))
  252. SmallestElement = Magnitude;
  253. ElementCount++;
  254. }
  255. /* Case where element is structurally zero */
  256. else
  257. { if (Data)
  258. printf(" ...");
  259. else
  260. putchar('.');
  261. }
  262. }
  263. putchar('\n');
  264. #if spCOMPLEX
  265. if (Matrix->Complex AND Data)
  266. { if (Header)
  267. printf(" ");
  268. for (J = StartCol; J <= StopCol; J++)
  269. { if (pImagElements[J - StartCol] != NULL)
  270. { printf(" %8.2gj",
  271. (double)pImagElements[J-StartCol]->Imag);
  272. }
  273. else printf(" ");
  274. }
  275. putchar('\n');
  276. }
  277. #endif /* spCOMPLEX */
  278. }
  279. /* Calculate index of first column in next group. */
  280. StartCol = StopCol;
  281. StartCol++;
  282. putchar('\n');
  283. }
  284. if (Header)
  285. { printf("\nLargest element in matrix = %-1.4g.\n", LargestElement);
  286. printf("Smallest element in matrix = %-1.4g.\n", SmallestElement);
  287. /* Search for largest and smallest diagonal values */
  288. for (I = 1; I <= Size; I++)
  289. { if (Matrix->Diag[I] != NULL)
  290. { Magnitude = ELEMENT_MAG( Matrix->Diag[I] );
  291. if ( Magnitude > LargestDiag ) LargestDiag = Magnitude;
  292. if ( Magnitude < SmallestDiag ) SmallestDiag = Magnitude;
  293. }
  294. }
  295. /* Print the largest and smallest diagonal values */
  296. if ( Matrix->Factored )
  297. { printf("\nLargest diagonal element = %-1.4g.\n", LargestDiag);
  298. printf("Smallest diagonal element = %-1.4g.\n", SmallestDiag);
  299. }
  300. else
  301. { printf("\nLargest pivot element = %-1.4g.\n", LargestDiag);
  302. printf("Smallest pivot element = %-1.4g.\n", SmallestDiag);
  303. }
  304. /* Calculate and print sparsity and number of fill-ins created. */
  305. printf("\nDensity = %2.2f%%.\n", ((double)ElementCount * 100.0)
  306. / (((double)Size * (double)Size)));
  307. if (NOT Matrix->NeedsOrdering)
  308. printf("Number of fill-ins = %1d.\n", Matrix->Fillins);
  309. }
  310. putchar('\n');
  311. (void)fflush(stdout);
  312. FREE(PrintOrdToIntColMap);
  313. FREE(PrintOrdToIntRowMap);
  314. return;
  315. }
  316. #endif /* DOCUMENTATION */
  317. /* Added to export the row and column maps to convert the
  318. internal matrix to an external form - JLM */
  319. void
  320. spRowColOrder(
  321. spMatrix eMatrix,
  322. int* OrdToIntRowMap,
  323. int* OrdToIntColMap
  324. )
  325. {
  326. MatrixPtr Matrix = (MatrixPtr)eMatrix;
  327. int I, Size;
  328. ASSERT_IS_SPARSE( Matrix );
  329. Size = Matrix->Size;
  330. if ( OrdToIntRowMap == NULL OR OrdToIntColMap == NULL)
  331. { Matrix->Error = spNO_MEMORY;
  332. return;
  333. }
  334. for (I = 1; I <= Size; I++)
  335. { OrdToIntRowMap[ Matrix->IntToExtRowMap[I] ] = I;
  336. OrdToIntColMap[ Matrix->IntToExtColMap[I] ] = I;
  337. }
  338. }