Radiation therapy is a clinical field in which computer simulations play a crucial role. Before patients undergo radiation therapy, an individual treatment plan for each patient needs to be created based on the specifics of their case (a process often referred to as treatment planning). The main aspect of this is setting the control parameters for the treatment machine so that the radioactive dose delivered to the patient is as concentrated to the tumor volume as possible. The inverse problem of determining such appropriate control parameters is typically formulated as an optimization problem, which, considering the complexity of modern treatment machines, requires computerized algorithms to solve accurately. Solving this optimization problem can be a key computational bottleneck in the treatment planning workflow. In many cases, finding a suitable treatment plan is a trial-and-error process, requiring multiple solutions of the optimization problems with different weights and parameters.

This thesis proposes different methods to enable the use of high performance computing (HPC) hardware to accelerate the optimization process in radiation therapy. We deal with two main computational aspects of the optimization workflow, the calculation of dose, gradients and objective functions; as well as the optimization solver itself. For dose calculation during optimization, we propose a CUDA kernel for sparse matrix-vector products tuned to dose deposition matrices from proton therapy. For the evaluation of the objective function -- which is often constructed as a weighted sum -- we develop a method to distribute the calculation of the objective function and its gradient across computational nodes using message passing.

For the optimization solver itself we first propose a task-based parallel implementation for Cholesky factorization of banded matrices. This can be an important kernel in interior point methods (IPM) for optimization, depending on the structure of the optimization constraints. The final two papers in this thesis deal with the adoption of iterative linear algebra in IPMs to enable GPU acceleration of the optimization solve. We develop an IPM using a doubly augmented formulation of the KKT linear system and Jacobi preconditioned conjugate gradient method, which we show can solve problems to acceptable accuracy as well as benefit substantially from GPU acceleration. Benchmarks against the commercial treatment planning system RayStation shows that our solver can improve performance by up to 4.4 times on some realistic cases.

Strålterapi är en klinisk gren där datorsimuleringar spelar en viktig roll. Innan patienters strålbehandlingar påbörjas behövs individuella behandlingsplaner skapas baserat på varje patientfalls specifika omständigheter (en process som ofta kallas dosplanering (eng. treament planning)). Huvudaspekten av denna process är att bestämma styrparametrar för behandlingsmaskinen så att stråldosen som levereras till patienten är så koncentrerad till tumören som möjligt. Det inversa problemet att bestämma sådana styrparametrar formuleras typiskt som ett optimeringsproblem, vilket, givet hur komplexa moderna behandlingsmaskiner är, kräver datoralgoritmer för att lösas precist. Att lösa detta optimeringsproblem kan vara en viktig flaskhals i dosplaneringsprocessen. I många fall är det en iterativ process att hitta en lämplig behandlingsplan som kräver att man testar och löser optimeringsproblemet med flera olika vikter och parametrar.

I denna avhandling studerar vi metoder för att möjliggöra användandet av hårdvara från högprestandaberäkningar för att accelerera optimeringsprocessen i strålterapi. Vi angriper två olika beräkningsaspekter i lösningen av optimeringsproblemet: beräkningen av stråldos, gradienter och målfunktioner, samt lösningen av optimeringsproblemet i sig. För dosberäkningen som krävs under optimeringen utvecklar vi en beräkningsmetod i CUDA för matris-vektor produkter med glesa matriser som är specifikt anpassad för dosmatriser (eng. dose deposition matrices) från protonterapi. För beräkningen av målfunktionen --- som i många fall är formulerad som en viktad summa --- utvecklar vi en metod för att distribuera beräkningen av målfunktionen och dess gradient över flera beräkningsnoder genom message passing.

För optimeringslösaren själv så utvecklar vi först en parallel implementation av Cholesky-faktorisering för bandade matriser med parallelisering baserad på tasks. Just den beräkningsalgoritmen kan vara viktig för många inrepunktsmetoder för optimering, baserat på strukturen av bivillkoren i optimeringsproblemet. De sista två publikationerna i avhandlingen studerar användandet av iterativa linjära lösare i inrepunktsmetoder för att möjliggöra GPU-accelerering av av optimeringslösaren. Vi utvecklar en inrepunksmetod med en dubbelaugmenterad formulering av KKT-systemet och en Jacobi-förkonditionerad CG-lösare (eng. conjugate gradient method). Vi demonstrerar att våran metod kan lösa optimeringsproblemen till tillräcklig noggranhet och dessutom ha stor nytta av GPU-acceleration. Prestandajämförelser med det kommersiella dosplaneringssystemet RayStation visar att våran lösare kan förbättra lösningstiden med upp till 4.4 gånger på några realistiska patientfall.