Speicherung von banded Matrix in Fortran

Question

Ich schrieb einen abgeleiteten Datentyp zum Speichern von banded Matrizen im komprimierten diagonalen Speicherformat; Insbesondere speichere ich jede Diagonale der Bandmatrix in einer Spalte der 2D-Arrays cds(1:N,-L:U), wo N die Anzahl der Zeilen des Vollmatrix und L und U ist die Anzahl der unteren und oberer Diagonalen (this question umfasst die Definition von der Typ).

Ich schrieb auch eine Funktion, um das Produkt zwischen einer Matrix in diesem CDS-Format und einem vollständigen Vektor durchzuführen. Um jedes Element des Produktvektors zu erhalten, werden die Elemente der entsprechenden Reihe von cds verwendet, die im Speicher nicht zusammenhängend sind, da die Sprache Fortran ist. Aus diesem Grund wanderte ich, wenn eine bessere Lösung wäre, die Diagonalen in der Reihe eines cds2(-L:U,1:N) Reihen zu speichern, die mir ziemlich vernünftig scheint.

Im Gegenteil here las ich

wir für die Matrix A ein Array val(1:n,-p:q) zuordnen kann. Die Deklaration mit umgekehrten Dimensionen (-p:q,n) entspricht dem LINPACK-Bandformat [132], das im Gegensatz zu komprimiertem Diagonalspeicher (CDS) keine effizient vektorisierbare Matrix-Vektor-Multiplikation zulässt, wenn p + q klein ist.

Was meiner Meinung nach gerade C angemessen erscheint. Was vermisse ich?

EDIT

Der Kern der Routine Performing Matrixvektorprodukte ist die folgende

DO i = A%lb(1), A%ub(1) 
    CDS_mat_x_full_vec(i) = DOT_PRODUCT(A%matrix(i,max(-lband,lv-i):min(uband,uv-i)), & 
             & v(max(i-lband,lv):min(i+uband,uv))) 
END DO 

(Wo lv und uv verwendet werden berücksichtigt den Fall des Vektors aus einem Index indiziert anders als 1.) Auf die Matrix A wird dann durch Zeilen zugegriffen.

Answer 1

Ich implementierte den abgeleiteten Typ, der die Diagonalen in einem Array val(-p:q,1:n) speichert und es ist schneller, wie ich vermutete. Also denke ich, dass der Link, auf den ich verwiesen habe, sich auf eine Speichersprache in einer Zeile als C und nicht auf eine Spalte als Fortran bezieht. (Oder es implementiert das Matrixprodukt in einer Weise, die ich mir nicht einmal vorstellen kann.)