Matrices gigantes

Siento haber tardado tanto en responder, pero es que respondí ayer noche y escribí un montón de cosas, ya lo creo. Pero he aquí mi sorpresa cuando veo que la pregunta sigue activa. Lo peor es que no me acuerdo de todas las soluciones que te propuse. Intentaré, resumidamente, reconstruir lo que pensé:
Al principio de todo te comentaba que, gracias a Dios te habías equivocado en los cálculos y eran en realidad 368 MB o algo así, y si tienes en cuenta que sólo hace falta almacenar un poco más de la mitad, por ser simétrica, quedaría en 181 MB. Además, en estos cálculos he tenido en cuenta que el tamaño de un double son 8 bytes y no 4.
Básicamente te propuse cinco soluciones:
1º Si tu matriz es esparcida (poco densa), es decir, con un porcentaje elevado de ceros, puedes almacenar sólo los valores no nulos junto con sus índices en una tabla de acceso rápico (tabla hash o árboles). Además para el cálculo de la matriz inversa podrías utilizar métodos iterativos muy eficientes para este tipo especial de matrices, los llamados métodos iterativos de solución de sistemas lineales de ecuaciones. Básicamente hay dos: método de Jacobi y método de Gauss-Seidel.
2º Es poco probable que tu matriz sea esparcida, por lo que tendrías que, para calcular la inversa de tu matriz, utilizar un método exacto de reducción de Gauss-Jordan o variaciones para reducir propagación de errores (en tu caso el uso de estas variaciones serían totalmente necesarias, al ser tu matriz gigantesca; estas técnicas se llaman de pivoteo). Se trataría de almacenar sólo las filas que se estén utilizando en los cálculos.
3º Manejar la forma diagonalizada de la matriz. Al ser simétrica, existe una matriz semejante que es diagonal, es decir, sólo necesitarías almacenar 10000 datos. Lo peor es que necesitarías almacenar 10000 autovectores con 10000 componentes cada uno. Esta opción no la he analizado muy bien. Tal vez sea la mejor porque puedes almacenar sólo 10000 ángulos de rotación y generar las matrices asemejantes en tiempo real, es decir, en total sólo necesitarías almacenar 20000 datos, pero se necesitarían muchos más cálculos. De todos modos, los tiempos de un disco duro son del orden de milisegundos, mientras que los de la memoria de nanosegundos y al final te saldría muy rentable los cálculos adicionales de generación sobre la marcha de autovectores.
4º Utilizar algún método de intercambio de bloques con caché en RAM. En este caso establecer una política óptima no sería muy difícil ya que sería específica para tu programa que actuaría de una forma conocida de antemano o a priori.
5ºSi puedes disponer de otros ordenadores conectados en red, puedes realizar cálculos paralelamente, además de que la cantidad de memoria total se incrementaría. En este caso, tendrías que utilizar algún tipo de protocolo (MPI) para comunicación de procesos distribuidos.
No dudes en preguntarme para más detalles de alguna de ellas. Mejor aún, si dispones de alguna otra información acerca de la matriz o cual es el propósito de la misma, puede ser útil para dar otras soluciones.

Hola Albertod
En primer lugar agradezco tu atención. Tu respuesta es completamente contundente e informada, se advierte con claridad tu dominio tanto en el área matemática como informática. He leído atentamente tus planteamientos y me llama mucho la atención los puntos 2 y 3 ya que como suponías mi matriz no es esparcida. Lo de la reducción Gauss-Jordan y las variaciones no lo conocía y suena bien ya que si se puedo reducir el tamaño de la matriz por alguno de estos métodos para invertirla sería fabuloso, ¿Me puedes informar sobre este tema?
Con respecto al punto 3 te solicito que si has averiguado algo porfavor me lo hagas saber, pero si no estonces solo cuentame lo del punto 2, la idea es no molestarte tanto con este tema porque imagino que debes tener otras preocupaciones más importantes. De todas maneras insisto en agradecer tus nobles y altruistas intenciones. Ahora te informo en más detalles lo que precisamente va a hacer mi programa (eso espero);
Tengo unas funciones escalares para los elemntos j <> j F(a,x,y,z)
y tengo una función escalar para los i=j W(a,x,y,z)
Bien con eso no hay problema ya que son las funciones que generan la matriz R,
Ahora el orden de la matriz R va a depender de la cantidad de elementos finitos en que desee dividir "el sistema de elctrodos" y que debe alcanzar como máximo 10.000x10.000 Ahí comienza el problema.
Ahora para obtener la "resistencia del sistema (r) " que es el objetivo del programa, debo invertir la matriz R y luego en esa matriz G (La inversa de R) sumar todas las filas, es decir, obtener un vector V de orden 1x10.000 donde el elemento v(1,1) de ese vector va a coresponder a g(1,1)+g(2,1)+g(3,1)+... g(n, 1) y el elemento v(1,2) va a coresponder a g(1,2)+g(2,2)+g(3,2)+... g(n, 2) y así.
Luego que obtengo tal vector V solo tengo que sumar todos sus elementos
s = v(1,1)+v(1,2)+v(1,3)+.......+v(1,n)
En otras palabras y en resumen, tengo que sumar todos los elementos de la matriz G
Y finalmente
r = 1/s
se ve sencillo no?
de una forma más elegante; r = [ [1]t x [R]-1 x [1] ]-1
Con [1]t = [1 1 1 1 ....... 1] de orden 1xn (la t es de traspuesta y el -1 de inverso)
Como puedes ver el problema principal consiste en invertir R. Espero haber expuesto con claridad los detalles del problema y esperaré tu respuesta con respecto a Gauss-Jordan y las variaciones. Muchas Gracias
Sinceramente
Rodrigo

El método de Gauss-Jordan es simplemente un método clásico de resolución de sistemas de ecuaciones. En el caso del cálculo de la matriz inversa, se parte de la matriz ampliada que es la matriz R junto con la matriz identidad a la derecha: [ R | I ] = A. Se trata de llegar a la configuración [ I | G], mediante pasos elementales que conducen a sistemas equivalentes. Al final G = R-1. En cada paso que desde 1 hasta n (n es la dimensión de R) se elije una fila i que tenga el elemento A(i, k) <> 0 (si no se encuentra este i, es que la matriz no tiene inversa). Se intercambian las filas i y k de A. Se divide la fila k por A(k, k). Ahora se producen ceros en la columna k, excepto en la fila k. Para eso, para cada i desde 1 hasta n distinto de k, se resta a la fila i el resultado de multiplicar A(i, k) por la fila k.
Las técnicas de pivoteo permiten reducir errores de redondeo que resultarían inadmisibles en algunos sistemas de ecuaciones. Para ello, la técnica más simple dice que en vez de buscar un i que tenga A(i, k) <> 0, buscamos un i que además tenga A(i, k) máximo. Esto servirá en la mayoría de los casos.
El hecho de que la matriz simétrica lo puedes aprovechar porque la inversa de una matriz simétrica es a su vez simétrica. No sé qué se puede sacar de esto exactamente para optimiar el algoritmo de Gauss, pero puede merecer la pena pensarlo.
Lo que no veo solución es en que tendrás que almacenar la matriz en disco. Cualquier simplificación parte de la matriz entera. Otra cosa es que la matriz sea siempre la misma para todas las ejecuciones del programa, entonces puedes utilizar técnicas de precondicionamiento de algoritmos.