Позвольте уточнить чуток, чтоб дать совет поточнее.
Правильно ли я понимаю, что у вас эрмитова разреженная матрица, вышедшая из дискретизации чего-то похожего на Хартри-Фока? В чем у вас матрица дана, есть ли предобуславливатель, каков он, на каком компьютере Вы это считаете?
Давидсон с предобуславливателем - самое оно, но тут как раз предобуславливатель - самое ключевое в этом гиблом деле звено. Причем правильнее тащить все подпространство во время решения, но K-ая итерация начинает поедать K*N слов памяти, где N - число неизвестных, поэтому используют ньютоновский предобуславливатель поверх обычного, чтоб уменьшить размер пространства.
По моему опыту на 80М в двойной точности скалярное произведение еще не сильно валится, особенно если делать реортогонализацию подпространства по два раза, но я для таких задач считаю скалярные произведения уже в 4-ой точности, а потом округляю до обычной двойной.
Коды старые на фортране есть, могу поделиться, что найду. Писал для своего диссера в 1999. Но код в общем-то сильно параллельный (на ту старую МВС100 или МВС1000 я уже запамятовал). Не факт, что все гладко полетит у Вас, да и потом мой проедобуславливатель на специальные конечные элементы был заточен, и там на 100М неизвестных давидсон за сотню итераций для 20-ой электронной пары сходился

На кривом базисе такие результаты получить не реально, и давидсон часто перестартовать надо будет.
Короче, расскажите про задачу чуток побольше, уверен, что смогу помочь.
Ой... Написал, а потом задумался... У Вас матрица 9К*9К? Не... влоб, переводя в плотную и далее лапаком в один вызов, и без вариантов. То что написано выше, это для матриц 80М*80М.
> У меня компилятор фортрана (и гфортран и пиджиай) ругаются, а точнее просто вылетают в ошибку, если размер массива превышает 81М (т.е. больше чем 9кХ9к).
а это - банально в фортране нельзя столько в стек засовывать, лечится тем, что память аллоцируется в Сишнике с помощью malloc, и передается в фортран.